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ABSTRACT 


This  thesis  explores  the  application  of  the  maximum  likelihood  parameter 
identification  technique  to  determine  the  wheel/rail  creep  coefficients  using 
dynamic  response  data.  The  equations  of  motion  for  the  dynamically  scaled 
wheelset  are  presented  and  the  reduced  form  of  the  maximum  likelihood  equa¬ 
tions  as  applicable  to  the  dynamically  scaled  wheelset  model  are  developed. 

The  maximum  likelihood  equations  were  formulated  into  a  maximum  likelihood 
algorithm  which  was  implemented  in  Fortran  IV.  Using  simulated  wheelset  data, 
the  effects  of  a  random  input  representation  of  the  track  versus  a  determin¬ 
istic  input  with  uncertainty  representation  are  determined.  The  effects  of 
various  levels  of  measurement  noise  are  also  examined.  This  preliminary 
analysis  indicates  that  the  deterministic  representation  of  the  track  input 
yields  better  results.  Representing  the  track  as  a  random  track  input  requires 
further  investigation  into  the  effects  of  longer  data  records  and  smaller  time 


steps  on  the  performance  of  the  maximum  likelihood  algorithm. 
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INTRODUCTION 


The  dynamic  response  of  railroad  vehicles  to  track  irregularity  inputs 
is  determined  by  interactions  between  the  vehicle  body  and  suspension  com¬ 
ponents,  kinematic  constraints  dependent  on  the  profiles  of  the  wheels  and 
rails,  and  friction  forces  generated  at  the  wheel/rail  interface.  These  lat¬ 
ter  forces,  known  as  creep  forces,  result  from  small  differential  velocities 
between  wheel  and  rail  in  directions  tangential  and  lateral  to  the  wheel. 
Linearized  dynamic  models  for  rail  vehicles  represent  creep  forces  using  the 
slope  of  the  creep  force  versus  relative  velocity  function  at  the  origin;  the 
slopes  are  known  in  the  literature  as  creep  coefficients. 

The  values  of  these  coefficients  has  been  determined  from  Hertzian  con¬ 
tact  theory  by  Kalker  (4)  and  by  measurements  under  steady-state  conditions 
(9.).  However,  previous  analysis  (12)  of  dynamic  response  data  has  indicated 
that  better  agreement  between  theory  and  experiment  is  obtained  if  values  for 
creep  coefficients  are  used  that  are  up  to  50 %  lower  than  determined  in  (4_) 
and  (9.).  The  "best"  values  for  creep  coefficients  under  dynamic  conditions 
have  not  been  determined  with  precision;  such  determination  would  be  of  great 
value  to  the  rail  research  technical  community. 

The  objective  of  this  thesis  is  the  estimation  of  creep  coefficients  from 
experimental  vehicle  response  data  using  advanced  parameter  identification 
techniques.  To  focus  the  effort  on  estimation  of  creep  coefficients,  data 
used  are  measurements  of  the  response  of  a  simplified  scale  model  vehicle, 


whose  characteristics  and  parameters  are  well-understood,  exclusive  of  the 
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creep  coefficients.  Use  of  the  simplified  model,  in  this  case  a  wheelset 
(single  axle  with  two  wheels  and  lateral  and  yaw  suspension,  capable  of  two 
degrees-of-freedom  motion),  permits  examination  of  the  fundamental  kinematic 
hunting  and  instability  characteristic  of  rail  vehicles,  without  the  effects 
of  more  complex  vehicle  configurations  obscuring  the  results.  The  use  of 
dynamically  scaled  models  (9_)  allows  better  control  of  experimental  condi¬ 
tions  than  expected  in  full  scale. 

The  method  used  to  estimate  the  creep  coefficients  is  the  maximum  likeli¬ 
hood  parameter  identification  technique.  This  technique  has  been  used  to 
identify  aircraft  stability  derivatives  successfully  for  several  years,  and 
represents  the  state-of-the-art  in  parameter  estimation.  The  specific  objec¬ 
tive  of  this  thesis  is  to  develop  and  implement  the  maximum  likelihood  method 
for  dynamically  scaled  wheelset  models.  In  addition  some  insight  into  the 
best  way  to  conduct  dynamic  experiments  involving  the  dynamically  scaled 
wheelset  model  was  desired. 

This  research  program  was  divided  into  two  phases.  The  first  was  the 
experimental  phase  during  which  the  dynamic  wheelset  experiments  were  con¬ 
ducted.  The  purpose  of  the  experiments  was  to  obtain  data  that  could  be 
processed  using  a  maximum  likelihood  processor.  The  second  stage  of  the 
research  project  was  a  developmental  one.  During  this  part  of  the  research 
program  the  maximum  likelihood  algorithm  was  developed  and  implemented  into 
a  Fortran  IV  computer  program.  The  final  phase  of  research  was  the  test  phase 
during  which  simulated  wheelset  data  was  generated  using  a  computer  model  of 
the  dynamically  scaled  wheelset  model.  The  maximum  likelihood  computer  pro¬ 
gram  was  then  checked  out  using  this  simulated  data. 

The  experimental  phase  was  conducted  first  so  that  as  much  information 
as  possible  about  the  wheelset  model  could  be  obtained  prior  to  the  development 
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of  the  maximum  likelihood  program.  The  actual  vheelset  data  taken  was  not 
analyzed  using  the  maximum  likelihood  program.  The  results  for  the  simu¬ 
lated  data  suggest  that  more  research  needs  to  be  done  with  simulated  data 
so  that  more  conclusive  results  can  be  drawn  when  the  actual  wheelset  data  is 
processed. 

The  body  of  this  thesis  is  divided  into  three  chapters.  The  first  pre¬ 
sents  the  generalized  maximum  likelihood  equations  and  then  develops  the 
reduced  equations  that  apply  to  the  dynamically  scaled  wheelset  model.  The 
second  chapter  discusses  the  experimental  phase  of  the  research  program. 

The  knowledge  gained  from  these  initial  experiments  will  be  used  to  develop 
an  improved  experimental  program  for  the  next  series  of  data.  The  final 
chapter  presents  the  results  obtained  for  the  simulated  data,  and  the  inter¬ 
pretation  of  these  results  in  terms  of  actual  wheelset  data. 


-***—  ay 


Chapter  1 


MAXIMUM  LIKELIHOOD  THEORY 


1.1  General  Form  of  the  Maximum  Likelihood  Equations 

The  maximum  likelihood  parameter  identification  technique  is  used  to 
identify  unknown  parameters  of  a  dynamic  system  given  measurements  of  some 
or  all  of  the  system's  state  variables.  This  method  maximizes  the  probability 
of  the  given  measurements  conditioned  on  the  vector  of  unknowns.  The  follow¬ 
ing  development  of  the  maximum  likelihood  equations  will  be  for  a  general 
case,  linear,  time-varying  system.  Except  where  otherwise  noted  the  equations 
in  this  section  were  taken  from  References  (8)  and  (3). 

The  state-space  form  of  the  differential  equations  for  a  linear  dynamic 
system  is 

x(t^  =  F(t)x_(t)  +  L(t)u(t)  +  G(t)w(t)  (l.l) 

where  x(t)  =  state  variable  vector 
u(t)  =  known  input  vector 
w(t)  =  input  disturbance  vector 
The  measurement  process  is  described  by  the  equation 

cXt)  =  H(t)x(t)  +  v(t)  (1.2) 

where  z_(t)  =  measurement  vector 

x(t)  =  state  variable  vector 
v(t)  =  measurement  noise  vector 

The  foundation  of  the  maximum  likelihood  method  is  maximizing  the  con¬ 
ditional  probability  of  the  vector  for  unknowns  given  the  measurements  of  the 


states 


5 


maximize  Plil*  ) 


where  p  denotes  the  probability  density 

0  =  vector  of  unknowns  to  be  identified 
zn  =  sequence  of  measurements 

Maximizing  p(0|zn)  is  equivalent  to  maximizing  the  log  likelihood  function. 

L  =  ln[p(8jzn)]  (1.3) 

where  L  A  log  likelihood  function 

Sequential  application  of  Baye's  Rule  to  the  conditional  probability  of  0  given 

z  yields: 
n 

p(ilzn)  =  p(znl£)  p(£lzn_i)  d.M 

=  p(znl£)  p(zn_1 1£)  p(£I  zn_2^ 

=  p(zn|0)  p(zn_x |_0)  p(zn_2 1 e.)  p(ilzn_3) 

N 

=  n  p(z.  | j9)  p(s_|o)  d.5) 

i=l  1 


Because  the  state  vector  is  dependent  upon  the  vector  of  unknowns  the 
equation 


N 

p(0jzn)  =  II  p(zi |£)  p(  0J 0 ) 


can  be  re-written  as 

p(£lzn) 


N 

n  p(z  |x.(0))  p(0jO) 
i=l  "”1 


(1.6) 


Assuming  the  input  disturbance  (w(t))  and  the  measurement  noise  (v(t)) 
are  gaussian,  then  the  formula  for  multivariate  gaussian  conditional  distri¬ 
bution  can  be  applied  because  p(z.  |x,(0_))  is  gaussian. 


p(z. |x-  (£) )  “ 


-E[z.])J 


b71(z. 
i  —a 


-  IfeiD 


(2ir) 


N/2 


lBil 


1/2 


(1.7) 
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In  equation  1.7  E[z^]  =  H_^  where  ^  is  the  output  of  a  Kalman  Filter.  The 

matrix  B.  is  defined  as  follows: 

1 


Bi  =  E[(zi  -  EIZJH.Z.  -  EfzJ)  ] 


B.  =  E[(z.  -  H.£. )(z.  -  H.fi. )  ] 

l  — l  l— 1  —i  ’  — 


B.  =  H.  P.  HT  +  R. 

l  ill  l 


(1 


Combining  equations  1.3,  1.6  and  1.7  yields 
N 


L  *  ""‘‘i"!  2„lN/2)  ,B  (1/2)  • 


-±(z.  -  H . x .  )  B.  (z. 

2  —l  l—i  l  —l 


H.S.  ) 

l—i 


]  •  p(9jo)}  (1 


L  =  -  4  Z  { ( z .  -  H.x. )T  B71(z.  -  H.fc. )  +  ln|B. | } 

2  .  .  -a  l-i  l  -l  l-a  '  i  1 

i=l 

-  (|)ln(27r)  +  ln[p(6|0)]  (1 

,  N  m  i 

L  =  -  —  Z  {(z.  -  H.S.  )  B.(z.  -H.5c.)+ln|B.|}  +  constant  (l 

2  ^  -l  1-1  i  — i  1-1  '  i  1 

The  maximum  likelihood  problem  involves  maximizing  equation  1.11  with 
respect  to  theta,  subject  to  the  constraints  of  the  Kalman  Filter  equations. 
The  Kalman  Filter  equations  must  be  satisfied  because  the  state  estimate  x^ 
and  the  state  covariance  (P^)  are  outputs  of  the  Kalman  Filter.  A  summary 
of  the  continuous  time,  time  varying  Kalman  Filter  equations  are  given  below 

C5). 

System  Model: 

x(t)  =  F(t)x(t)  +  L(t)u(t)  +  G(t)wCt)  (l 

Measurement  Model: 

z_(t)  =  H(t)x(t)  +  v(t)  (l 

State  Estimates: 

kit)  =  F(t)i(t)  +  L(t)uCt)  +  K(t)  z{t)  -  H(t)$(t)  (1 


.8) 

.9) 

.10) 

.11) 


.12) 

.13) 

.lM 
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Error  Covariance  Propagation: 

P(t)  =  FCt)P(t)  +  P(t)FT(t)  +  G(t)Q(t)GT(t)  -  P(t)HT(t)R_1(t)H(t)P(t)  (1.15) 

Kalman  Gain  Matrix: 

K(t)  =  P(t)HT(t)R-1(t)  (1.16) 

where  v(t)  N(0,R(t))  (1.17) 

and  w(t)  'V'  N(0,Q(t))  (l.l8) 


The  notation  used  in  equation  1.17  denotes  that  v(t)  is  a  random  se¬ 
quence  with  a  gaussian  probability  density  with  a  mean  value  and  variance 
defined  below. 

E[v(t)]  =  0  (1.19) 

E[(_n+.)  -  E[v(t )  ] )  (v(t)  -  E[v(t)])T]  =  R(t)  (1.20) 

Equation  1. .  •,**  is  solved  on  the  digital  computer  for  discrete  values  of 
using  equati  •  1.21  through  1.2**. 

xtt  +At)  -  $(t  ,At)2(t  )  +  0(t  )z(t  )  +  T(t  )u(t  )  ( 1 . 21 ) 

—  O  O—O  O—O  O—O 

where  [F(t  )  -  K(t  )H(t  )]At  4>(t,At)  =  e  (1.22) 

o  o  o 

0(t  )  =  <J>(t  ,At)[F(t  )  -  K(t  )H(t  )]_1  •  [I  -  $_1(t  , At ) ] K ( t  )  (1.23) 

o  o  o  o  O  0  0 

r(t  )  =  4>(t  ,At)[F(t  )  -  K(t  )H(t  )]_1  •  [I  -  4>-1(t  ,At)]L(t  )  (1.2**) 

o  o  o  o  o  o  o 

The  error  covariance  propagation  equation  (1.15)  is  solved  using  the 

Kalraan-Engler  solution. 


P(t  +At )  =  [<)>,  (t  ,At)  +  <t)iT1(t  ,At)P(t  )] 
o  Ay  o  aa  o  o 


[^VAt)  +  ^X(to,At)P(to)] 


-1 


where 


jr(tQ+At) 

A(t  +  t) 

U“  o  J 

V 


>1 

• 

)] 

A(tol 

(1.25) 


(1.26) 


T 


8 


<f>  (t  .At)  <p  At  ,At) 

*(t  ,At )  =  ™  0  yX  0 

Ky(VAt)  *AX(VAt)- 

The  equations  for  £_(t)  and  A_(t)  are  given  in  equation  1.28. 


(1.27) 


-FT(t) 


X(t ) I  G(t)Q(t)G  (t) 


HT(t)R-1(t)H(t)l  f  jr(t 


(1.28) 


Therefore 


^(t^At)  =  exp 


-FT(t  ) 
o 

G(t  )Q(t.)GT(t  ) 
o  o  o 


HT(t  )R-1(t  )H(t  ) 

O  0  0 

F(t  ) 
o 


(1.29) 


The  linear  system  of  equations  in  1.29  is  derived  in  Reference  ( 3_) . 

The  likelihood  function  can  be  expanded  into  a  Taylor  series  with  only 
three  terms  because  the  likelihood  function  for  Gaussian  conditional  proba¬ 
bility  distribution  is  quadratic. 


L(0)  =  L(0  )  + 

o  90 


(0-0  )  +  (0-0  )T 

i  a^L 

0  o-o 

2  982 

0 

— 0 

_ 

-o 

(1.30) 


Taking  the  partial  derivative  of  L(9_)  with  respect  to  0_  and  setting  it 
equal  to  zero  to  find  the  maximum: 


=  3L 

9*  9- 

e 

*4 

90^ 

(9*-8  )  =  0 

—  — o 


(1.31) 


where  is  the  value  of  8_  that  results  in  the  maximum  of  the  likelihood  func¬ 
tion.  Solving  for  9_*  yields: 

l\2r ,  l-i  rT  i  IT 


0«  =  9  - 
—  — o 


362  9 


(1.32) 


In  equation  1.32  A0,  or  the  step  in  theta  is  equal  to: 


9 


A6  =  - 


302  6 
-o. 


8ie 

— o 


(1.33) 


■gg’  is  defined,  as  the  gradient  and  is  equal  to: 

31  N  3v(t  ) 

—  i=l  — 


where 


3B(t  ) 

~2V(ti)B  ~W~B  1(ti)v(ti) 


.  ,  3B(t.) 

*  2  tr  B'  <V  ~lf- 


v(t.)  =  z(t.)  -  H ( t ^ ) x_( t ^ ) 


(1.3U) 


(1.35) 


3v(t . )  3£(t.)  3H(t . ) 

4  =  TJ  (  4-  )  ^  ^ 

36  38  -  36 


(1.36) 


(V  r 

5T"  “  H(ti>  [p(ti 


3(H  (t  ))  9P(t  ) 

)  - i -  +  - L_ 

1  30  30 


V  T 

jT"H 


3H(t. )  m  3R(t. ) 

+  -3T-p(ti)H  (V  +  30 


(1.37) 


The  term  — r  is  defined  as  the  Fisher  information  matrix. 
3(T 

„2.  H  3vT(t.)  ,  3v(t . ) 


m  .  9B(t  )  3v(t  ) 

-  (2)v  (ql®  <V  -jj-B-htj)  -Jg3- 


-  itr 


B-1 ( t  ± 


9B(t,  ) 

>  -ST"  B' 


3B(t. )' 

)  - 1 — 

i  36 


3fi(t,) 


(1.38) 


The  term  — ^ in  equation  1.36  is  referred  to  as  the  sensitivity  term 

and  is  derived  from  the  Kalman  Filter  state  estimate  equation  (1.1*0. 
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3x(t ) 


36 


-  F(t) 


afcU)  3F(t)  aL(t)  3i(t) 

-^0—  +  i(t)  +  u(t)  -  K(t)H(t)  - 


30 


-mi  ait) i(t)  .an £(t) _ 


3K(t) 


30  T  30  -Vl"  30 

See  Appendix  A  for  the  derivation  of  equation  1.39. 

Equation  1.39  can  be  solved  for  discrete  values  of  — gg —  on  the  digital 
computer  using  the  following  equations 


H(t)*(t) 

3x(t ) 


(1.39) 


where 


3x(t^.  )  3S_(t  +At)  3£(t_) 

A  - S -  =  4>(t  ,  t) 


30 


°  +  Alt  )*(t  ) 


30  ”"'0’  *"  30 


+  'i'lt  )u(t  )  +  T(t  )z(t  ) 
o  o  o  —  o 

[Fit  )-K(t  )H(t  ) ]At 
*(t  ,At )  =  e  °  °  ° 


o—o 


(1.1*0) 

(1.1*1) 


A(tQ)  =  #(to,At)[F(tQ)  -  K(to)H(to)] 


-1 


[I  -  <D_1(tQ,At)] 


3F(t  ) 
o 

30 


-  «V 


3H(t  ) 
o 

3K(t  ) 
o 

30 

'  30 

K(t  )H(t  )1  1. 

o 

0 

3L(t  )1 
0 

3£ 

Hit  ) 
o 


[I  -  $_1(t  ,At)] 
o 


T(t  )  =  <*>(t  , At )  [F(t  )  -  Kit  )H(t  )] 
o  o  o  o  o 


-1 


(1.1*2) 


(1.1*3) 


[I  -  4>"1(to,  t)] 


3K(t  ) 
o 

30 


3pCti ) 

The  term  — jjg —  in  equation  1.37  is  solved  using  equation  1.1*5. 


3P(t  +  t) 
o 


30 


■[ 


4>  (t  ,At)  +  <j>  (t  ,At) 
yy  o’  yy  o’ 


3P(t _) 


r  3PU  ) 

4’yy(to’At)  +  V/VAt)  30 


30 

-l 


VI 

>0 


(1.1*1*) 


(1.1*5) 


9P(t  +  t) 

Equation  1.1*5  is  an  iterative  solution  to  - -  and  is  derived  in 

00 

Appendix  B. 

The  final  unknown  term  in  the  sensitivity  equation  is  ^  —  .  This  term 
is  derived  by  differentiating  equation  1.16  with  respect  to  theta 


9K(t) 

90 


=  P(t) 


9HT(t) 

90 


R-1 (t )  +  HT(t) 


90 


+ 


9P(t ) 
90 


HT(t)R-1(t) 


(1.1*6) 


1.2  Application  of  the  Maximum  Likelihood  Equations  to  the  Wheelset 


The  equations  developed  in  the  last  section  apply  to  any  linear  dynamic 
system.  In  this  section  those  equations  will  be  tailored  to  the  dynamically 
scaled  wheelset  problem. 

The  wheelset  equations  of  motion  are  given  in  equations  1.1*7  through  1.50. 


where 


x(t)  =  Fx(t )  +  Gw(t ! 


F  = 


.Of  +  f.y 
A22  Ct 

mv 


'22 


K  7 


2V' 


av 

5? 


G  = 


x  = 


(y-6 ) 


w  =  6 


(1.1*7) 


(1.1*6) 


0 

L-U 

and  w(t)  *v  N(0,Q) 

See  Reference  (ll_)  for  a  detailed  development  of  the  wheelset  equations  of 
motion.  Appendix  C  contains  definitions  for  the  variable  of  the  F  and  G 


(1.1*9) 


(1.50) 


matrix. 
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The  wheelset  measurement  equation  is  given  in  1.51 

z_(  t )  =  H  x_(  t )  +  v  v  t ) 

where  100 

H  =  F  0  1  ol 


U.51/ 


0  0  1  -1 

and  v(t)  a-  !i(0,R)  U.53) 

The  system  dynamics  matrix  (F),  the  system  random  input  matrix  (0),  and  the 
system  measurement  matrix  (K),  the  covariance  matrix  for  the  random  input 
disturbance  (Q)  and  the  covariance  matrix  for  the  measurement  noise  (R)  are 
all  time-invariant. 

The  creep  coefficients  f  and  f^  in  the  system  dynamics  matrix  are  the 
parameters  to  be  identified  by  the  maximum  likelihood  method. 

For  the  purposes  of  simplifying  the  maxim-urn  likelihood  equations  the 
actual  creep  coefficients  will  be  defined  as  a  multiple  of  their  nominal 
values . 

p  e  O  f  :  1  c  1 

*00  ■  OO  'OO  —  •  > 


r*  _ 

‘22 

S22 

f 

22 

S'  _ 

6  „ 

*11 

A11 

lx 

s  „  = 

6 

=  s 

22 

11 

where,  for  simplicity  S0^  =  8^  =  6  (l.5o' 

The  nominal  values  for  the  creep  coefficients  have  been  identified  in  other 
research  on  railroad  vehicle  stability  (M. 

Combining  equations  1.1*6  and  1.51*  through  i.5o  the  system  dynamics  matrix 
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Since  the  system  dynamics  matrix  contains  the  only  parameter  to  be 
identified  the  vector  of  unknowns  9_  reduces  to  a  scalar. 

6=6  (1.55) 

The  maximum  likelihood  equations  of  the  last  sect-ion  can  be  greatly 
simplified  for  two  reasons: 

1)  the  wheelset  is  a  time-invariant  system. 

2)  the  system  dynamics  matrix  (F)  contains  the  only  parameter  to  be  identified. 
A  summary  of  the  time-invariant  Kalman  Filter  Equations  for  the  wheelset 


system  is  given  below 
State  Estimates: 

x(t)  =  Fx(t)  +  Kz[ t )  -  KHx(t)  (1.59) 

Error  Covariance  Propagation: 

P(t)  =  FP(t )  +  P(t)FT  +  9QGT  -  P(t)  H'“'R~1HP(t )  (1.60 

Kalman  Gain  Matrix: 

K  =  P  hV1  (1.61) 

ss 

V(t)  't  rUO.R)  (1.62) 

w(t)  ~  H(0,Q)  (1.63) 

Equation  1.21,  the  solution  to  the  Kalman  filter  estimates  equation 
(l.lM  ,  reduces  to 

x(t  +At)  =  $(At)  x(t  )  +  0  z(t  )  (1.62) 

—  o  —  o  —  o 

where  . , . [F-KKlAt  ,, 

$(At)  =  e L  1  (l.o5) 

0  =  4>(At)  [F-KH]-1  •  [I  -  1  ( At ) ] K  (1.66) 

Pss  in  equation  1.6l  is  the  steady  state  solution  to  equation  1.6? 

P(tQ+At)  =  [<i^(At)  +  <|>^(At)  P(tQ)]  *  [4>yy(At)  +  4>^(At)  P(to)]_1  ( 1.6T ) 


This  equation  is  the  same  as  equation  1.25  except  that  is  no  longer  dependent 


Ik 


on  t  ,  but  is  a  function  only  of  At,  as  shown  in  equation  1.68. 


4>(At)  =  exp 


-F 


m  -« 

HXR  H 


At 


(1.68) 


^GQG  F 

Because  all  of  the  elements  of  the  matrix  in  equation  1.68  are  time  invariant, 
the  iterative  solution  for  P  has  to  be  solved  onlv  once.  This  value  of  P 

S'j 

is  used  to  calculate  the  Kalman  gain  matrix  according  to  equation  l.Cl.  The 
Kalman  gain  matrix  for  a  time  invariant  system  such  as  the  wheelset  remains 
constant  for  the  entire  Kalman  gain  problem. 

Taking  into  account  the  two  simplifications  listed  above,  the  maximum 
likelihood  equations  reduce  to: 


L  =  -  -  [  Z  v'vt.)B  v It.)  ♦  In |B j ]  (] 

i=l 

The  constant  term  in  equation  l.i:  be  dropped  because  it  shifts  the 

likelihood  function  up  or  down  but  does  not  change  the  location  of  the  maxi¬ 
mum  of  the  likelihood  function. 

IJ 


,69) 


|  =  ^  ^  -  f  V^ti)B-1  |  B-1v(t1)  +  i  trtB-1  f  ] 


9v( t . ) 

l 

39 


3x ( t .  ) 

i 


39 


3B  _  3P  1 

39  -  39 


3  L 
392 


K  3vT(t.)B-1  3v(t. 


=  £ 
i=l 


90—-  B_1 


39 


3v(t .  ) 

i 

39 


(1.T0) 

(1.71) 

(1.72) 


Maximum  Likelihood  Algorithm 


The  equations  of  the  preceding  section  were  implemented  into  a  Fortran 
IV  computer  program.  This  section  discusses  the  algorithm  used  to  solve  the 
wheelset  maximum  likelihood  problem  on  the  digital  computer. 

The  maximum  likelihood  method  is  a  batch  processor  in  that  the  entire 
set  of  measurements  is  used  for  each  iteration  of  the  maximum  likelihood 
equations  (8_) .  Each  iteration  of  the  maximum  likelihood  algorithm  produces 
an  estimate  of  the  value  of  the  parameter  to  be  identified  which  will  maximize 
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the  likelihood  function  given  in  equation  I.69.  This  estimate  0_*  then 
replaces  0^  within  the  program  and  another  iteration  is  made.  When  the  fol¬ 
lowing  condition  is  met,  iteration  is  stopped. 

6*  -  0  =  A6  <  e  (l.8l) 

—  ^0  — 

During  the  development  phase  of  a  maximum  likelihood  processor  it  is 

helpful  to  calculate  the  value  of  the  likelihood  function  (L),  the  gradient 

2 

Dij  3  L 

(t-j-)  ,  and  the  Fisher  information  matrix  ( — — )  for  a  range  of  6  where  the 

3i  302 

maximum  of  the  likelihood  function  is  expected  to  lie.  Changes  in  the  flow¬ 
chart  of  Figure  1.1  are  shown  in  Figure  1.2.  The  flowchart  at  the  bottom  of 
the  dotted  box  replaces  the  A0  test  at  the  end  of  the  flowchart  in  Figure  1.1. 
This  procedure  is  beneficial  in  that  the  maximum  likelihood  urogram  can  be 
checked  out  without  the  risk  of  an  infinite  loop  being  set  up  because  of  a 
convergence  problem. 

Ho  mention  has  been  made  so  far  of  the  method  used  to  calculate  the  state 
transition  matrix  (STM).  On  a  digital  computer,  a  very  fast  and  accurate  way 
to  calculate  the  STM  is  through  a  power  series  expansion. 

$(  At )  =  I  +  AAt  +  ^-A2(At)2  +  jj-A3(At)3  +  ...  +  A  (1.82) 

The  series  is  truncated  when  the  following  condition  is  satisfied 
N+l  ...  II  .  .  -| 

(  I  A1  ( At ) 1  —+  I)  -  (  l  A1  (At)1  rj  +  I)  <  e  (1.33) 

i=l  lp  i=l  1- 

where  e  is  some  error  criterion  chosen  by  the  user.  A  value  of 

e  -  1.0  •  10-10  (1.8’-) 

is  usually  considered  to  be  sufficient.  A  listing  of  the  maximum  likelihood 
program  is  contained  in  Appendix  E. 


WWHi  1 


Figure  1.1  Maximum  Likelihood  Algorithm  (8.). 
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1.1  Test  Data  Generation 

In  order  to  ensure  that  the  maximum  likelihood  computer  program  imple¬ 
ment  the  maximum  likelihood  algorithm  properly,  it  was  tested  using  simulated 
wheelset  data.  The  data  was  generated  using  the  state  transition  matrix  ap¬ 
proach  to  solve  the  following  equation. 

x  =  Fx  +  Gw  (1.85) 

where  F  and  G  are  specified  by  equations  1.57  and  I.I9  respectively.  Due  to 
limitations  in  computer  storage  a  total  of  1000  data  points  for  each  state  was 
all  that  could  be  generated. 

The  time  step  (At)  was  chosen  by  examining  the  eigenvalues  of  the  F  matrix. 
Typical  values  for  the  poles  of  the  wheelset  model  are: 

a1  =  -210.7  sec  1 

a2  3  =  -0.U23  ±  J  5.982  sec"1 

These  poles  indicate  that  transient  responses  are  characterized  by  a  lightly 

damped  sinusoid  (oon  =  6.0  sec  1;  C,  =  0.07)  plus  a  rapidly  converging  exponential 

decay  (t  =  =  1.7  msec).  In  the  frequency  domain  a  sharp  resonant  peak  occurs 

°1 

at  the  natural  frequency,  with  an  additional  breakpoint  at  u  8  1/t. 

As  the  estimates  of  the  creep  coefficients  change,  the  damping  ratio  of 
the  low  frequency  resonance  changes,  as  does  the  location  of  the  higher  fre¬ 
quency  real  pole.  In  selecting  the  time  interval  AT  and  the  record  length  NAT 
it  is  desirable  for  information  on  both  high  and  low  frequency  modes  to  be 
present  in  the  data.  To  obtain  six  data  points  per  cycle  at  the  higher  fre¬ 
quency  f  =  o^/2tt  a  AT  is  chosen  to  be  0.005  sec.  Since  computer  storage 
constraints  limited  the  record  length  to  N  =  1000  points,  the  record  length  is 
5  seconds,  which  corresponds  to  1.8  cycles  of  the  low  frequency  mode. 
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Large  values  for  At  were  tried  but  the  state  transition  matrix  approach  used 
in  the  maximum  likelihood  program  had  convergence  problems.  A  At  =  .005  was 
the  largest  time  step  that  could  be  used  without  causing  numerical  problems. 

In  addition  to  the  random  track  input  w  in  equation  1.85,  there  was  also 
measurement  noise  added  to  the  state  vectors. 

z_  =  Hx_  +  v  (1.86) 

w  and  v  were  generated  by  summing  10  random  number  vectors  which  were  genera¬ 
ted  by  an  APL*  operator.  The  random  vectors  generated  by  successive  runs 
of  this  APL  operator  are  not  correlated.  Ten  of  these  random  vectors  were 
summed  to  force  the  distribution  of  u  and  w  to  be  approximately  Gaussian  (1_) . 

Once  these  gaussian  random  vectors  were  formed  they  could  be  manipulated  to 
force  their  means  and  mean  squared  values  to  equal  what  the  user  desired. 

Although  maximum  likelihood  theory  requires  w  and  v  to  be  white  noise  be¬ 
cause  these  vectors  were  of  finite  length  they  were  bandwidth  limited.  Also 
since  only  10  random  vectors  were  summed,  the  probability  distributions  of  w 
and  v  were  not  exactly  gaussian,  only  approximately  gaussian.  A  listing  of 
the  APL  functions  used  to  generate  the  simulated  wheelset  data  is  contained  in 
Appendix  F. 


*  APL:  A  Programming  Language.  Implemented  on  the  Princeton  University  IBM 
370/158  time-sharing  system. 
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Chapter  2 


EXPERIMENTAL  PROGRAM 


2.1  0b,1  ective 

The  experimental  program  was  designed  to  make  measurements  of  the  wheel- 
set  state  variables.  The  third  order  wheelset' model  has  the  following  state 
variables : 

v  =  lateral  velocity 
=  yaw  angle 

( y— 6 )  =  wheelset  to  rail  centerline  relative  lateral  displacement 
See  Reference  (ll)  for  a  detailed  development  of  the  wheelset  equations  of 
motion. 

2.2  Experimental  Facility 

The  entire  experimental  program  was  carried  out  on  the  Princeton  Dynamic 
Model  Track.  Previous  research  at  the  Dynamic  Model  Track  involved  the 
development  and  validation  of  the  dynamically  scaled  wheelset  model  and  the 
measurement,  restraint  and  propulsion  system  for  the  wheelset  (9.). 

The  experimental  aoparatus  consisted  of  the  wheelset  model  and  the  idler 
carriage.  Both  of  these  components  rode  on  the  kOO  foot  LEXAI.’  track.  The 
idler  carriage  surrounded  the  wheelset  but  was  dynamically  isolated  from  it. 
The  idler  carriage  provided  restraints  to  prevent  the  wheelset  from  completely 
derailing.  The  idler  carriage  was  also  part  of  the  linkage  that  connects  the 
wheelset  to  the  propulsion  unit.  Figures  2.1  and  2.2  diagram  the  wheelset- 
idler  carriage  relationship. 
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The  propulsion  system  was  a  hydraulically  operated  drive  unit  that  rode 
along  a  steel  I-beam  guideway  above  the  LEXAIJ  track.  The  hydraulic  drive 
unit  had  a  feedback  control  system  for  maintaining  a  constant,  user-selected 
velocity. 

Suspended  from  the  propulsion  unit  on  instrument  racks  were  the  on-board 
analog  computer,  the  bridge  amplifiers,  and  the  portable  four  track  FM  analog 
cassette  recorder.  The  bridge  amplifiers  were  used  to  amplify  the  output  of 
the  lateral  accelerometer. 

2.3  Transducers 

Figure  2.3  is  a  diagram  of  transducer  placement  on  the  wheelset  and  Figure 
2.1  traces  the  output  of  the  transducers  through  the  simal  conditioninr 
equipment  to  the  recorder. 

Lateral  velocity  could  not  be  measured  directly,  therefore  an  accelero¬ 
meter  was  used  to  measure  lateral  acceleration.  A  bridge  amplifier  was  used 
to  excite  the  differential  type  accelerometer  and  to  amplify  the  output  of 
the  accelerometer.  The  lateral  acceleration  signal  will  have  to  be  integrated 
on  a  digital  computer  to  obtain  lateral  velocity. 

The  yaw  angle  was  measured  directly  using  a  potentiometer.  The  equation 
relating  potentiometer  output  (v)  to  the  vaw  angle  (v)  is: 

0  =  mv  +  b  (2.1) 

where  m  is  the  slope  of  the  straight  line  described  by  this  function  and  b  is 
the  y-intercept.  The  intercept  b  varied  from  dav  to  day  due  to  amplifier  drift. 
A  pre-run  and  post-run  procedure  was  developed  so  that  the  value  of  b  could  be 
determined  at  the  start  of  every  day. 
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The  variations  in  yaw  angle  during  a  run  were  smaller  in  magnitude  than 
the  value  of  the  y-intercept.  In  order  to  obtain  the  best  resolution  on  the 
cassette  tape  an  electrical  bias  was  summed  with  the  yaw  potentiometer  output 
on  the  on-board  analog  computer.  The  value  of  the  electrical  bias  was  ad¬ 
justed  from  day  to  day  so  that  the  output  of  the  analog  computer  circuit 
equalled  zero  when  the  yaw  angle  was  zero.  The  equation  that  describes  the 
analog  computer  output  is: 

H  =  ^  +  e  (2.2) 

m 

where  q  =  the  output  of  the  analog  computer 
e  =  the  electrical  bias 

A  displacement  transducer  (DCDT)  was  mounted  on  each  side  of  the  wheelset. 
Each  displacement  transducer  measured  the  distance  of  the  wheelset  from  the 
inside  edge  of  the  track.  Figure  2.5  diagrams  the  relationship  between  the 
variable  measured  by  the  DCDT's  and  the  state  variable  (y-5).  The  equations 
for  the  DCDT's  are: 

V  =  m  (v-6  )  (2.3) 

r  t  d 

=  ra^  (y-5-^)  (2.M 

where  =  voltage  output  of  the  right  DCDT 

VT  =  voltage  output  of  the  left  DCDT 

u 

m^  =  slope  of  the  straight  line  function 

m  =  slope  of  the  straight  line  function 

Li 

rn 

Multiplying  VT  by  —  and  adding  V  to  VT : 

m_  r  ii 

L 

m 

V  +  — 
r  rr^ 

This  is  the  signal  that  was  recorded  during  the  experiment.  The  on-board 


V  =  m  (y— 6  )  +  m  (y-6  ) 
L  r  r  r  L 


(2.5) 
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analog  computer  was  used  to  multiply  the  left  DCDT  signal  by  —  .  The  state 
_  m 

variable  (y-6)  can  be  derived  from  the  recorded  signal  bv  dividing  V  +  —  VT 

r  “l  L 

by  2mr  . 


V  +  —  VT  ,  -  r  . 

r  L  2  y  -  mr  (6  +  6  ) 


(2.6) 


(y-6)  =  y  - 


(6l+62} 

=  y  -  6 


(2.7) 


The  DCDT's  have  a  linear  range  of  operation  which  is  much  smaller  than 
their  total  range  of  operation.  The  linear  range  of  operation  is  slightly 
different  for  each  DCDT  as  shown  in  Figures  2.6  and  2.7.  Each  DCDT  was  posi¬ 
tioned  so  that  its  output  was  the  center  of  its  linear  range  when  the  wheelset 
was  centered  on  the  track.  Because  the  DCDT's  were  not  mounted  symmetrically 
about  the  wheelset 's  longitudinal  centerline,  a  bias  was  introduced  into  the 
measuring  process. 


m^  (y-6j  )  +  b 


(2.2) 


V  =  m  (y-6  ) 
r  r  J  r 


(2.o) 


where  b  =  position  bias 

The  same  algebraic  operations  as  previously  yield: 


mm  m 

^VT  =  rn  (y-6.  )  + 

"l  Ij  "l  l  l  "l 


(2.10) 


V  =  m  (y-6  ) 
r  r  r 


(2.11) 


m  m 

V  +  —  V.  =  m  (y-6T )  +  —  b  +  m  (y-5  ) 
rn^L  r  J  L  rJr 


(2.12) 


=  2m  y  -  m  (6+6  )  +  —  b 
r  r  L  r 


(2.13) 
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In  addition  to  the  three  variables  mentioned  above,  the  wheelset  axle  angular 

velocity  U  =  v/r  )  in  the  rolling  direction  was  also  measured.  Although  the  angular 
o 

velocity  is  not  a  state  variable,  it  is  necessary  to  construct  the  velocity 
dependent  state  matrix  F.  The  output  of  the  tachometer  was  recorded  directly 
with  no  processing  on  the  on-board  analog  computer  as  with  the  lateral  displace¬ 
ment  and  yaw  angle  signals.  Sample  responses  from  the  above  transducers  are 
given  in  Figure  2.8. 

2.1*  Experimental  Procedure 

All  of  the  above  mentioned  signals  were  recorded  on  a  portable,  four 
track  FM  analog  cassette  type  recorder*  which  was  carried  on  board  the 
drive  unit.  The  recorder  was  started  and  stopped  manually  before  and  after 
each  run. 

The  bridge  amplifier  that  powered  the  accelerometer  was  balanced  at  the 
beginning  of  each  day  to  eliminate  the  drift  in  amplifier  output  that  occurred 
from  day  to  day. 

The  on-board  analog  computer  had  a  built-in  calibration  circuit  for  the 
accelerometer  bridge  amplifier.  This  calibration  circuit  produced  a  known 
voltage  that  appeared  to  the  accelerometer  as  an  acceleration  load.  This 
acceleration  load  remained  constant  from  day  to  day.  In  order  to  determine 
the  input-output  relationship  of  the  accelerometer,  the  accelerometer  output 


*  Philips  Mini  Log  U  Portable  Analog  Cassette  Recorder. 
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was  recorded  during  the  calibration.  This  calibration  test  was  performed 
before  every  data  run  to  account  for  amplifier  drift  between  data  runs. 

2.5  Measurement  Noise 

There  were  two  sources  of  uncertainty  in  the  data.  The  first  of  these 
was  electrical  noise.  The  bridge  amplifier,  analog  computer  and  hydraulic 
pump  for  the  propulsion  unit  all  operate  on  U00  cycle  alternating  current. 

Vibration  of  the  cassette  recorder  during  data  runs  was  another  source 
of  noisy  data.  The  steel  I-beam  guideway  that  supported  the  hydraulic  drive 
unit  had  gaps  between  the  I-beams.  The  vibration  that  resulted  when  the  drive 
unit  wheels  hit  these  gaps  was  transmitted  to  the  recorder,  and  showed  up  as 
noise  on  the  data  tapes.  As  mentioned  earlier,  the  recorder  was  suspended  from 
the  drive  unit  on  an  instrument  rack  t/hich  had  shock  absorbers  to  help  iso¬ 
late  the  recorder  from  vibration,  nonetheless,  some  vibrations  still  affected 
the  recorder  and  added  noise  to  the  data. 

In  addition  to  the  two  definite  sources  of  noise,  there  is  a  third  possible 
source  of  noise.  It  is  known  that  there  are  small  gaps  between  lengths  of 
LEXAD  rail.  These  small  gaps  are  the  result  of  contraction  and  expansion  of 
the  LEXAU  material  due  to  temperature  changes.  The  800  foot  building  which 
houses  the  Princeton  Dynamic  Model  Track  does  not  have  a  completely  controlled 
environment;  therefore,  temperature  variations  of  as  much  as  50°F  are  possible. 
Although  no  ruantitative  analysis  has  been  carried  out  to  verify  this  hypo¬ 
theses,  it  is  possible  for  the  wheelset  structural  modes  to  be  excited  when 
the  wheels  hit  these  gaps.  The  transducers  mounted  on  the  wheelset,  especially 
the  lateral  accelerometer,  could  detect  this  vibration,  thereby  adding  noise 


to  the  data 


Finally  to  insure  the  validity  of  the  data,  the  LEXAN  track  and  wheels 
were  cleaned  before  every  run.  Dirt  and  dust  on  the  track  or  wheels  could 
alter  rail/wheel  adhesion,  in  which  case  the  wheelset  model  would  be  invalid. 
Absolute  methanol  was  used  as  the  cleaning  agent  because  it  left  no  residue 
on  the  track. 

2.6  Filtering 

To  remove  the  noise  due  to  the  electrical  system  and  vibration  of  the 
wheelset  structural  modes  all  of  the  signals  were  filtered  off-line.  The 
acceleration  signal  was  filtered  with  a  second  order  low-pass  filter.  The 
yaw  angle  and  lateral  displacement  signals  were  filtered  using  seventh  or 
eighth  order  band-pass  filters,  and  the  rotational  velocity  signal  was  filter¬ 
ed  using  a  first  order  low-pass  filter.  The  cutoff  frequencies  for  the  fil¬ 
ters  are  given  in  Table  2.1. 

Table  2.1 

Filter  Cutoff  Frequencies 


Signal 

Filter  Type 

High  Pass 
Cutoff 

Low  Pass 
Cutoff 

Lateral  Acceleration 

low-pass 

N/A 

b9.9  Hz 

Yaw  Angle 

band-pass 

.1  Hz 

60  Hz 

Rolling  Angular  Velocity 

low-pass 

N/A 

20  Hz 

Wheelset  to  Rail  Relative 
Lateral  Displacement 

band-pass 

.1  Hz 

60  Hz 

2.7  Digitization 


After  filtering  the  data  had  to  be  digitized  so  it  could  be  stored  and 
analyzed  on  a  digital  computer.  As  mentioned  in  the  previous  section,  the 
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lowest  permissible  sampling  rate  was  determined  when  the  low-pass  filter 
cutoff  frequencies  were  chosen.  It  was  decided  to  use  a  sampling  rate  of 
120  samples /second  per  channel.  The  A/D  converter  used  for  digitization 
was  a  12  bit  sample  and  hold  A/D  converter.*  This  A/D  converter  was  controlled 
by  a  minicomputer  system.**  The  minicomputer  had  two  memories  for  storing 
data  (memory  A  and  memory  B)  each  of  which  had  the  capacity  to  store  512  data 
points.  When  the  digitizing  process  started,  the  computer  stored  the  digi¬ 
tized  data  in  memory  A.  Once  memory  A  was  full  the  computer  began  filling 
memory  B  while  at  the  same  time  writing  the  data  in  memory  A  to  magnetic  tape. 
Since  the  A/D  converter  was  a  sample  and  hold  type,  there  was  only  a  several 
nanosecond  time  delay  between  the  sample  for  each  channel.  However,  the  switch¬ 
ing  process  from  memory  A  to  memory  B  took  3  msec.,  therefore  the  sampling 
interval  between  the  512th  data  point  and  the  513th  data  point  was  11.33  msec. 
The  sampling  interval  between  all  other  sets  of  four  data  points  was  8.33 
msec.  Figure  2.8  depicts  this  shift  in  sampling  interval.  Appendix  G  is  a 
listing  of  the  program  that  controlled  the  A/D  converter  during  the  digitiz¬ 
ing  process. 


*  Preston  GMAD-1  Analog  to  Digital  Conversion  System. 

**  Hewlett  Packard  HP  1000  System. 


2  3  4 


0.00833  SEC. 
(NOMINAL  AT  FOR 
120  SAMPLES /SEC) 


5  *126  127  128  129  130  131  132  SAMPLE 

NUMBER 


O.OII33  SEC.  0  00833  SEC. 

(NOMINAL  AT  +  (NOMINAL  AT) 

0.003  SEC  SWITCH 
FROM  MEMORY  A 
TO  B) 
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Chapter  3 
RESULTS 


3.1  Outline  of  Test  Cases 

In  this  chapter  the  results  obtained  from  the  maximum  likelihood  computer 
program  using  the  simulated  data  will  be  presented.  As  mentioned  in  Chapter 
1,  the  AFL  program  which  generated  the  simulated  data  had  provisions  for  any 
combination  of  deterministic  input,  random  input,  initial  condition  and  mea¬ 
surement  noise.  Table  3.1  is  a  summary  of  the  maximum  livelihood  program 
testing. 

3.1.1  Explanation  of  Test  Case  Data 

The  LEXAiJ  track  lateral  alignment  is  ideally  a  white  r.cise  velocity  input 
with  a  gaussian  probability  density.  Because  of  these  characteristics  the 
track  input  does  not  have  to  be  measured  but  can  be  described  as  a  random 
input  into  the  wheelset  model.  In  this  case,  it  is  the  Q  matrix  which  gives 
the  maximum  likelihood  processor  all  of  its  information  about  the  track  input. 
It  is  also  possible  to  measure  the  track  input,  although  there  would  be  some 
uncertainty  involved.  If  the  track  input  were  measured,  then  the  measured 
values  for  the  track  input  would  become  a  deterministic  input  and  the  uncertai 
ty  in  the  track  measurements  would  be  described  as  a  rando:  track  input. 

For  the  first  seven  test  cases  the  track  input  is  treated  as  a  random 
track  input,  with  the  only  information  describing  it  contained  in  the  Q  matrix 

In  test  cases  7  and  8  the  track  input  is  treated  as  a  deterministic  input 
The  random  track  input  for  these  cases  represents  the  uncertainty  in  the 
deterministic  track  input. 
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In  terms  of  generating  the  simulated  vheelset  data,  the  two  represen¬ 
tations  of  the  track  input  are  handled  in  the  following  way.  For  the  case 
when  the  track  input  is  considered  to  be  a  totally  random  input,  a  random 
vector  is  used  as  the  input  to  the  syster  model  which  generates  the  data. 

When  this  data  is  processed,  the  maximum  likelihood  program  (specifically  the 
Kalman  Filter)  is  not  tol  i  what  the*  run  i-  vector  was,  the  only  information 
it  is  given  is  the  covariance  matrix  for  the  vector  Q. 

When  the  track  input  is  considered  to  be  deterministic,  the  system  model 
which  generates  the  data  has  two  random  vectors  as  inputs.  One  of  these  ran¬ 
dom  vectors  represents  the  deterministic  input; the  other  represents  the  uncer¬ 
tainty  in  the  deterministic  input.  This  time  when  the  data  generated  by  the 
system  model  are  analyzed  by  the  maximum  likelihood  processor,  the  Kalman  Filte 
is  given  the  vector  that  represents  the  deterministic  input,  but  it  is  given 
only  the  covariance  matrix  for  the  vector  that  represents  the  uncertainty  in 
the  deterministic  input.  Figure  3.1  diagrams  the  generation  of  the  simulated 
data  for  both  representations  of  the  track  input.  This  figure  is  not  designe 
to  represent  the  entire  process  for  generating  simulated  data  (4here  can  be 
measurement  noise)  nor  is  it  designed  to  represent  all  of  the  inputs  that  go 
into  the  maximum  likelihood  processor.  Its  only  purpose  is  to  demonstrate  the 
variations  in  representing  the  track  input. 

3.1.2  Purpose  of  Each  Test  Case 

As  stated  previously,  test  cases  8  and  9  are  different  from  the  other 
seven  test  cases  in  how  the  track  input  is  represented.  Test  case  9  is  dif¬ 
ferent  from  test  case  8  in  that  the  level  of  uncertainty  in  the  deterministic 
track  input  is  much  higher.  These  two  cases  give  some  indication  of  how  well 
the  real  track  input  would  have  to  be  measured  to  significantly  increase  the 


deterministic  track  input  with  some  uncertainty 


performance  of  the  maximum  likelihood  processor  over  the  case  when  the  track 


is  considered  to  be  a  random  track  input. 

Test  cases  1  through  5  were  run  to  find  out  how  various  combinations  of 
initial  condition,  random  track  input  and  measurement  noise  affect  the  maximum 
likelihood  processor.  Test  case  1  through  3  were  also  run  to  find  out  how 
well  that  maximum  likelihood  processor  performs  when  it  is  given  incorrect 
information  about  the  random  track  input  and  the  measurement  noise.  In  other 
words  the  simulated  data  for  test  cases  1  through  3  was  generated  with  a 
random  track  input  and  measurement  noise  that  had  statistics  as  given  in 
Table  3.1.  However,  when  this  data  was  analyzed  using  the  maximum  likelihood 
processor,  the  Kalman  Filter  was  given  incorrect  information  about  the  sta¬ 
tistics  of  the  w  and  v_  used  to  generate  the  data.  This  incorrect  information 
took  the  form  of  an  incorrect  Q  and  R  matrix.  Table  3.2  shows  the  true  Q  and 
R  matrices  for  w  and  v  versus  the  Q  and  R  matrices  actually  used  by  the  Kalman 
Filter . 


3.2  Method  of  Testing  the  Maximum  Likelihood  Parameter 


Essentially  the  maximum  likelihood  program  lias  3  main  parts.  The  first 
of  these  calculates  the  likelihood  function  based  on  the  equation 

li 


L  =  -  |  £  (v1(t.)B  1v(t.)  +  In  | B | ) 

i=l 


(3.1) 


The  likelihood  function  has  two  main  components  as  can  be  seen  in  equation  3.1. 

IJ 

L  ..  = -  ^  E  vT(t. )B_1v(t. )  (3.2) 

observation  2  .  .  1  1 

i=l 


The  second  part  of  the  maximum  likelihood  equation  calculates 

which  is  referred  to  as  the  gradient,  and  the  third  part  of  the  program  cal- 
o 
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culates  — —  known  as  the  Fisher  Information  Matrix.  The  inverse  of  the 

dQ_ 

Fisher  Information  Matrix  is  defined  as  the  Cramer-Rao  Lower  Bound  on  the 
covariance  of  the  0  estimates.  The  maximum  likelihood  approaches  the  lower 
bound  asymptotically  (8).  The  gradient  and  the  Fisher  Information  Matrix  are 
used  to  calculate  the  estimate  of  the  value  of  the  parameter  to  be  identified, 
according  to  the  equation 


where 


0*  =  0  -  A0 

—  — o 

ao  ..-Ir^T  J  „  92L 
A0  =  M  [  gQ-J  and  M  =  — ^ 

—  90 


As  stated  in  Chapter  1,  the  maximum  likelihood  technique  involves  iterating 

2 

8L  3  L 

through  equations  for  L,  rr  ,  — —  ,  and  0*  until  A0  becomes  smaller  than  some 

902 

error  criterion  (e).  According  to  maximum  likelihood  theory  as  e  -*■  0 ,  6*  ap¬ 
proaches  the  true  identity  of  the  unknown  parameter.  At  the  true  value  for 
0_*  the  likelihood  function  has  its  maximum,  and  according  to  maximum  likeli¬ 
hood  theory,  for  the  technique  to  converge  the  likelihood  function  must  have 
a  maximum.  Because  it  was  not  certain  that  the  likelihood  function  for  all 
the  test  cases  had  a  maximum,  the  maximum  likelihood  processor  was  not  allowed 

to  converge.  Instead  the  likelihood  function  was  given  specific  values  of 

2 

3L  3  l 

theta  for  which  to  calculate  L,  Cb^as>  and  — -g  .  Theta  ranged  iron  .50  to 

—  90 

1.45  in  increments  of  .05.  In  all  of  the  simulated  data  generated,  the  fol¬ 
lowing  value  for  the  unknown  parameter  was  used: 

B  =  1.00 

therefore  the  peak  in  the  likelihood  function  for  all  cases  should  be  at 


0  =  8  =  1.00 


For  each  value  of  theta  in  the  range  specified  above,  the  maximum  likelihood 


] 


program  iterated  through  its  equations  once  calculating: 


(1) 

L 

(2) 

f 

bias 

(3) 

^observation 

(M 

3L 

9£ 

_2 

(5) 

3  L 

2 

30* 

(6) 

A0 

(7) 

6* 

Figure  3.2  is  a  sample  of  the  computer  program  output  for  one  iteration. 
After  the  maximum  likelihood  program  finished  the  iteration  for  0  =  1.^5 
the  five  quantities  listed  above  are  plotted  for  the  range  of  theta. 


3.3  Test  Case  Results 

A  summary  of  the  resul+s  of  the  maximum  likelihood  program  for  test 
cases  1  through  9  are  given  in  Table  3.3.  These  results  will  be  discussed 
in  three  sections.  The  first  section  will  discuss  the  likelihood  function 
for  each  test  case.  The  second  section  will  discuss  the  results  for  the 
gradient  and  the  Fisher  Information  Matrix  and  the  third  will  examine  the 
application  of  the  test  data  results  to  the  analysis  of  real  wheelset  data. 

3.3.1  Likelihood  Function 

In  this  section  the  test  cases  will  be  discussed  in  three  groups  as  indi¬ 
cated  on  the  right  hand  side  of  Table  3.1. 


Examining  test  cases  8  and  9  first.  Table  3.1  shows  that  a  deterministic 
track  input  as  well  as  a  random  track  input  was  used  to  generate  the  simulated 
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Figure  3.2  Computer  output  for  the  last  iteration  of  test  case  number  8. 
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Table  3.3 

Summary  of  Results  for  Test  Cases 


Test  Case 

Likelihood  Function 
Maximum  6  = 

Observation  Term 
Maximum  0  = 

Bias  Term 
Maximum  0  = 

1 

.95 

.95 

N/A 

2 

N/A 

.95 

N/A 

3 

.95 

.95 

IJ/A 

h 

N/A 

N/A 

N/A 

5 

.95 

.95 

N/A 

6 

N/A 

N/A 

N/A 

7 

.65 

N/A 

N/A 

8 

.85 

1.10 

N/A 

9 

N/A 

N/A 

N/A 

NOTE:  N/A 

=  not  applicable,  no  maximum 

occurs  or  only  a  local 

maximum  occurs. 
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wheelset  data.  As  presented  in  Chapter  1,  the  wheelset  equations  of  notion 
in  state  vector  form  are: 

x.  =  Fx_  +  Gw  (3.5) 

where  w  =  random  track  input 

The  Kalman  Filter  state  estimates  for  this  case  are  given  by  the  equation 

x  =  (F-KH)x  +  Kz  (3.6) 


o  incorporate  deterministic  and  random  track  inputs  the  following  changes 
in  the  above  two  equations  are  necessary. 

x_  =  Fx.  +  Gu  +  Gw  (3.7) 

where  u  =  deterministic  track  input 

w  =  random  track  input  due  to  the  uncertainty  in  u 

x  =  (F-KH)x  +  Kz_  +  Gu  (3.8) 

Although  u  is  a  white  noise  input  it  is  called  deterministic  because  it  was 
an  input  into  the  system  model  which  generated  the  data  (equation  3.7)  and 
it  was  an  input  into  the  Kalman  Filter  of  the  maximum  likelihood  program 
(equation  3.3).  w  is  called  a  random  track  input  because  it  appears  as  an 
input  into  the  wheelset  system  (equation  3.7)  but  it  is  not  an  input  into  the 
Kalman  Filter  (equation  3.3).  The  Kalman  Filter  in  the  maximum  likelihood 
program  has  the  Q  matrix  as  its  only  source  of  information  about  w. 

The  difference  between  test  cases  8  and  9  is  the  size  of  the  random  track 
input  as  compared  to  the  deterministic  track  input.  Case  9  has  a  larger  ran¬ 
dom  track  input  than  case  8,  which  signifies  that  case  9  has  more  uncertainty 
in  the  deterministic  track  input  than  case  8.  As  Table  3.3  shows,  for  case  8 
the  likelihood  function  had  a  maximum  but  for  case  9  the  likelihood  function 
did  not  have  a  maximum.  This  would  suggest  that  the  maximum  likelihood  method 
needs  a  well-defined  deterministic  track  input  if  it  is  to  identify  a  peak  in 
the  likelihood  function.  However,  maximum  likelihood  theory  has  no  restrictions 


52 


This  suggests  that  the  maximum  likelihood  processor  may  need  more  measure¬ 
ments — a  longer  data  record — when  the  deterministic  track  input  has  a  large 
amount  of  uncertainty.  In  reference  (6_)  an  increase  in  the  number  of  obser¬ 
vations  is  shown  to  cause  a  more  pronounced  maximum  in  the  likelihood  function 
Although  there  are  no  quantitative  results  which  prove  more  measurements 
will  produce  a  peak  in  the  likelihood  function  for  the  dynamically  scaled 
wneelset  case,  this  is  an  area  in  which  further  research  should  be  done. 

For  all  of  the  test  cases  in  group  1,  the  Kalman  Filter  in  the  maximum 
likelihood  program  was  given  the  following  values  for  Q  and  E,  as  shown  ir. 
Table  3.1  and  Table  3.2. 


Q  =  1.0 


R  = 


.01 

0 

0 


0 

.01 

0 


0 

0 

.01 


These  five  test  cases  were  run  for  diagnostic  purposes  but  have  been  included 
here  because  they  exhibit  some  interesting  trends. 

Case  2  is  called  the  perfect  observation  case  in  the  literature.  In 
reference  (T_)  the  likelihood  function  is  shown  to  reduce  to  equation  3.0  for 
this  case. 


L 


xi  T  —1 

I  (t. )Q  1  (t. ) 
i  i 


(  3. 


where  v(t.)  =  z(t.)  -  0z(t.  )  (3. 

i  i  l-l 

when  Q  is  known.  Essentially  equation  3.9  says  that  the  likelihood  function 
simplifies  to  the  observation  term  for  the  perfect  measurement  case.  The 
maximum  likelihood  processor  did  not  analyze  this  limiting  case  because 
equations  3.9  and  3.10  were  not  incorporated  into  the  maximum  likelihood 
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algorithm  for  test  case  number  2.  The  results  of  test  case  2  do  provide  in¬ 
sight  into  the  performance  of  maximum  likelihood  processor  when  it  has  in¬ 
correct  or  imprecise  information  about  w  or  v_„  Figures  3.3,  3.^  and  3.5 
contain  plots  of  the  output  of  the  maximum  likelihood  program  for  test  case  2. 

Test  case  2  has  a  random  track  input  only  which  makes  it  a  limiting  case 
for  the  effect  of  ;u  versus  w  that  was  examined  in  cases  8  and  9.  Test  "'ase  2 
suggests  that  the  more  randomness  there  is  in  the  track  input,  the  more  dif¬ 
ficulty  the  maximum  likelihood  processor  has  in  identifying  maximum  for  the 
likelihood  function.  There  are  no  theoretical  limitations  of  this  type  f:r 
maximum  likelihood.  As  in  test  cases  8  and  9  one  practical  asnect  of  imrle- 
menting  the  maximum  likelihood  technique  that  should  be  investigated  is  data 
record  length.  Although  several  assumptions  were  made  in  the  generation  of 
the  simulated  data  as  discussed  in  section  l.k  it  is  believed  that  these 
assumptions  will  not  seriously  affect  the  performance  of  the  maximum  likeli¬ 
hood  processor. 

Case  number  3  is  a  special  case  of  maximum  likelihood  identification 
because  there  is  no  random  track  input.  The  likelihood  function  reduces  to 
the  observation  term  for  this  case  also  ;however  ;the  weighting  matrix  is  dif¬ 
ferent  than  in  the  no  measurement  noise  case. 

L  =  E  v^t.  )R-1v(t.  )  ( 3. 

i=l  1  1 

where  v(t.)  =  z(t. )  -  Hx(t.)  (3. 

ii  l 

when  R  is  known.  In  this  limiting  case  a  Kalman  Filter  is  used  by  the  maximum 
likelihood  processor  as  shown  in  equation  3.12,  however  the  weighting  matrix 
is  not  dependent  on  P  but  only  on  R.  The  observation  term  calculated  by  the 
maximum  likelihood  program  for  test  case  3  was 


Figure  3.1»  Plot  of  the  observation 


Figure  3.5  Plot  of  the  bian  term  for  t'  at  cane 
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1  ^  T  1 

K  =  -  =;  2  vi(t.)B'1v(t.) 

(3.13) 

i=l  1 

where 

B  =  HPH1  +  R 

(3.1M 

and 

vCt^ )  =  z(ti  )  -  Hx(ti ) 

(3.15) 

The  state  estimates  as  calculated  by  the  maximum  likelihood  program  are  not  as 
accurate  as  they  could  be  because  the  Kalman  Filter  was  given  an  incorrect 
value  for  Q  (Table  3.2).  Also  the  weighting  matrix  used  in  equation  3.13  is 
very  different  from  the  weighting  matrix  used  in  equation  3.11.  Because  the 
maximum  likelihood  algorithm  for  test  case  3  did  not  incorporate  equation  3.11 
and  because  Q  was  incorrect,  the  limiting  case  of  no  process  noise  was  not 
analyzed  correctly. 

Test  case  1  is  a  combination  of  the  two  limiting  eases  presented  above. 
The  likelihood  function  reduces  to: 


i 


necessarily  those  that  would  be  obtained  if  the  no  random  track  input-no 
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measurement  noise  case  were  implemented  correctly  using  equations  3.1o  and 
3.17. 

Test  case  U  and  test  case  5  are  the  only  members  of  group  1  for  which 
the  maximum  likelihood  program  was  given  correct  values  for  Q  and  R.  The 
results  of  test  cases  k  and  5  given.  As  demonstrated  in  other  test  cases  the 
maximum  likelihood  processor  did  not  identify  a  peak  in  the  likelihood  function 
when  the  only  input  driving  the  wheelset  system  is  a  random  track  input.  From 
the  test  cases  performed  as  part  of  this  research  program  no  conclusive  expla¬ 
nations  can  be  given  for  the  above  problem.  As  suggested  before,  a  longer 
data  record  in  these  cases  may  be  necessary  however  there  is  no  direct  evidence 
to  support  this  conclusion. 

Test  cases  6  and  7  were  run  to  see  if  the  ratio  of  Q  to  R  nai  any  affect 
on  the  performance  of  the  maximum  likelihood  processor.  For  test  case  6  the 
ratio  of  Q  to  R  is  the  same  as  that  for  test  cases  ^  and  5.  According  to 
maximum  likelihood  if  the  values  of  Q  and  R  are  changed  but  their  ratio  remains 
constant  then  the  observation  term  will  remain  the  same  but  the  bias  term 
will  change.  Graphs  of  the  observation  term  for  test  case  l  and  test  case  6 
are  presented  in  Figure  3.6  and  3.7  respectively.  Comparing  these  two  plots 
shows  that  the  observation  term  did  not  change  between  these  two  cases.  Com¬ 
paring  Figures  3.3  and  3.9  it  can  be  seen  that  the  bias  term  did  change  its 
range  of  values,  but  not  its  general  shape.  The  likelihood  functions  for 
case  and  case  6  are  given  in  Figures  3.10  and  3.11.  The  likelihood  function 
has  the  same  general  shape  for  each  case,  however,  the  location  of  the  local 
maximum  does  change.  This  agreement  between  test  cases  A  and  6  further  supports 
the  validity  of  the  maximum  likelihood  algorithm  given  in  Chapter  1,  and  the 
computer  program  which  implements  it. 
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igure  3.6  Plot  of  the  observation  term  for  test  case  h. 
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Figure  3.7  Plot  of  the  observation  term  for  tent  case 


M  J».  2  JO 


Figure  3.8  Flot  of  the  bias  term  for  test  case  h 


Figure  3.10  Plot  of  the  likelihood  function  for  test  cns 


Figure  3.11  Plot  of  the  likelihood  function  for  tent  cane 
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Test  case  6  has  as  its  only  input  a  randon  input  disturbance  and  in 
accordance  with  the  results  of  other  similar  test  cases,  its  likelihood 
function  does  not  have  a  maximum.  Test  case  7  also  has  a  random  disturbance 
as  the  only  input  to  the  wheelset  model;  however,  the  likelihood  function  in 
the  case  7  results  does  have  a  maximum  at  0  =  .65.  A  plot  of  this  likeli¬ 
hood  function  is  given  in  Figure  3.12.  The  results  of  test  case  7  indicate 
that  the  test  data  record  was  long  enough  for  the  maximum  likelihood  processor 
to  identify  a  peak  in  the  likelihood  function,  although  the  peak  occurred  at 
0_  =  .65  instead  of  _0_  =  1.00.  There  is  no  indication  that  a  longer  record 
will  shift  the  peak  from  9  =  .65  to  9  =  1.00,  although  this  is  an  area  which 
needs  further  research.  The  fact  that  the  measurement  noise  in  case  7  is 
much  smaller  than  in  case  6  may  be  the  reason  the  maximum  likelihood  processcr 
was  able  to  identify  a  peak  in  the  likelihood  function. 

There  is  another  possible  explanation  for  the  problems  the  likelihood 
function  had  in  identifying  a  peak  in  the  likelihood  function.  As  discussed 
in  section  l.U  a  At  =  .005  sec  was  used  to  generate  the  test  data.  This  time 
step  resulted  in  6  data  points  per  cycle  of  the  highest  frequency  response  of 
the  wheelset.  The  maximum  likelihood  processor  may  need  more  information 
about  this  particular  mode  of  response  of  the  wheelset.  This  could  be  accom¬ 
plished  by  using  a  smaller  At.  However,  because  the  nresent  method  used  to 
generate  the  simulated  data  is  restricted  to  1000  data  points  per  state  varia¬ 
ble,  a  smaller  At  would  produce  data  with  fewer  cycles  of  the  low  frequency 
response  of  the  wheelset.  In  order  to  investigate  the  effects  of  a  smaller  At 
a  new  method  for  generating  the  simulated  data  will  have  to  be  developed. 

Since  all  of  the  simulated  wheelset  data  was  generated  with  the  same  At, 
and  for  several  cases  a  peak  in  the  likelihood  function  was  identified  there 
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kelihood  function  for  tent  ease 
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does  not  seem  to  be  a  need  for  a  smaller  At.  However,  in  those  cases  where 
the  simulated  data  was  generated  using  a  random  track  input,  the  maximum 
likelihood  processor  was  given  minimal  information  about  the  input  into  the 
wheelset  model  that  generated  the  data — only  the  covariance  matrix  for  the 
random  track  input.  To  make  up  for  this  limited  information  about  the  input 
the  maximum  likelihood  processor  may  need  better  information  in  other  areas. 
This  "better"  information  includes  more  data  points  per  cycle  of  the  high  fre¬ 
quency  response  mode  of  the  wheelset — obtained  by  a  smaller  At — and  more  ob¬ 
servations  of  the  low  frequency  response  of  the  wheelset — obtained  with  a 
larger  data  record. 
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n 


=  ^3  *  m32 
3  2 

mi  =  m21 

“U  =  mU3 


(3.20) 

(3.21) 

(3.22) 


Equations  3.18  through  3.22  were  used  to  calculate  the  derivatives  of  the 
likelihood  function  and  the  gradient.  The  slope  at  the  first  data  point  was 
set  equal  to  the  slope  between  the  first  two  data  points.  The  slope  at  the 
last  data  point  was  set  equal  to  the  slope  between  the  last  two  data  points. 
This  method  of  averaging  slopes  works  very  well  when  the  data  points  form  a 
smooth  function.  In  many  cases  the  likelihood  function  data  points  and  the 
gradient  data  points  do  not  form  smooth  functions  and  the  derivatives  calcu¬ 
lated  using  this  second  method  can  have  large  errors.  Since  the  only  objective 
was  to  check  the  reasonableness  of  the  gradient  and  second  partial  data,  the 
second  method  of  averaging  slopes  was  used.  In  general  the  difference  between 
the  gradient  and  the  derivative  of  L  as  determined  by  the  average  slope  method 
is  on  the  order  of  5/«.  The  difference  between  the  second  partial  and  the 
slope  of  the  gradient  as  determined  by  the  averaging  method  was  considerably 
larger — on  the  order  of  35u.  This  difference  for  the  second  partial  was  con¬ 
siderably  larger  for  test  cases  4  and  6.  The  reason  the  difference  between  the 
second  partial  and  the  slope  of  the  gradient  may  be  so  large  for  test  cases  4 
an i  o  is  that  determining  derivatives  by  averaging  slopes  produces  large  errors 
■  iu":  points  which  do  not  form  smooth  functions.  The  data  points  for  the 
•  •'  -  -aser.  “  and  6  are  not  as  smooth  as  in  other  cases.  See  Figures 
.  a  re  iiscrepar.cv  is  most  probably  due  to  errors  in  the 

.  Information  Matrix  rather  than  in  the  calculation 

•  ma  i  lent  and  the  second  partial  was 
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not  carried  out  for  test  cases  6  and  9  because  based  on  the  results  for  test 
cases  1  through  7  it  was  concluded  that  the  quasilinearization  section  of  the 
maximum  likelihood  algorithm  was  implemented  correctly  in  the  maximum  likeli¬ 
hood  computer  program. 

3.  A  Application  of  Results  to  Actual  Wheelset  Data 

The  test  data  cases  were  analyzed  by  the  maximum  likelihood  program  for 
the  following  reasons: 

(1)  to  ensure  that  the  computer  program  was  working  properly. 

(2)  to  identify  any  potential  problems  that  may  arise  when  the  actual  wheelse 
data  is  analyzed. 

(3)  to  make  recommendations  for  further  research,  based  on  the  conclusions 
reached  in  (2)  above. 

The  first  objective  was  attained  as  documented  in  the  previous  sections.  The 
purpose  cf  this  section  is  to  discuss  items  an:  '3). 

Ir.  every  case  except  one  where  the  only  input  in*'  the  wheelset  motel 
was  the  ran ior.  track  input  the  likelihood  function  did  n  t  have  a  maximum, 
r^y^rvcil  Airies  in  the  or  evi  u.j  3  6  ?  1 1  r  r.  r'ir.t  r.r,.o  o",  -.-ert-f*  i  *  o  t  *.  o  like!  i  — 
h  oo  i  fun.-t  i;r,  failed  tc  have  a  maximum  becuase  ‘he  data  re’rr :  nee  Je  :  *~  bv 
longer.  The  test  cases,  run  as  part  cf  this  research  program  dr  not  form  a 
:arge  enough  experimental  base  to  state  the  above  print  with  any  certainty. 

It  is  suggested  that  further  research  in  this  area  be  directed  to  find  out  how 
the  length  of  the  data  record  affects  the  performance  of  the  maximum  likeli¬ 
hood  processor.  As  mentioned  in  Chapter  1,  the  present  method  of  generating 
simulated  data  is  limited  to  1000  observations  because  of  limitations  in 
computer  storage,  therefore  another  method  of  generating  simulated  data  will 


have  to  be  designed. 

Another  area  of  research  which  should  be  examined  is  the  effect  of  a 
smaller  At  in  generating  the  test  data.  A  smaller  At  would  provide  more 
information  about  the  high  frequency  response  mode  of  the  wheelset.  In  order 
not  to  lose  information  about  the  low  frequency  mode,  data  records  longer 
than  1000  points  per  state  variable  will  have  to  be  used.  As  pointed  out 
above  this  requires  a  new  method  of  generating  simulated  data. 

Appendix  H  contains  plots  of  the  likelihood  function,  observation  time 
and  bias  tern  for  test  cases  1  through  9. 


As  stated  in  the  Introduction,  the  objective  of  this  research  program  was 
the  development  and  implementation  of  a  maximum  likelihood  parameter  identi¬ 
fication  algorithm  applicable  to  a  dynamically  scaled  wheelset  model.  The 
following  is  a  summary  of  the  conclusions  that  can  be  drawn  from  the  results 
of  this  research  program. 

l)  The  maximum  likelihood  parameter  identification  equations  can  be  tailored 
to  the  problem  of  identifying  creep  coefficients  for  the  dynamically 
scaled  wheelset  model. 

T)  The  reasonable  results  obtained  for  many  of  the  test  cases  proves  that 
the  implementation  of  the  maximum  likelihood  algrrithr.  in  Fortran  1'.' 
was  performed  correctly. 

:■)  More  research  involving  simulated  wheelset  data  is  necessary  bef're 
actual  wheelset  data  can  be  processed.  It  is  recommended  that  the 
following  areas  be  investigated. 

a  >  the  effect  of  data  record  length  on  the  performance  of  the  maxi  mum 
likelihood  processor 

b )  advantages  in  representing  the  track  input  as  a  deterministic 
input  with  uncertainty  as  compared  to  representing  the  track 
input  as  a  random  track  input 

c)  the  effect  of  a  smaller  At  on  the  performance  of  the  maximum 
likelihood  processor. 

The  primary  motivation  for  the  above  three  recommendations  is  that  it 


seems  reasonable  to  assume  the  maximum  likelihood  processor  works  better 
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vhen  it  has  more  information  about  the  system  and  about  the  inputs  into  the 
system.  If  the  track  input  is  a  random  track  input  then  the  covariance 
matrix  is  the  only  information  the  maximum  likelihood  processor  receives 
about  this  random  track  input.  In  order  to  make  up  for  this  limited  informa¬ 
tion  about  the  input  into  the  wheelset  system,  the  maximum  likelihood  pro¬ 
cessor  should  receive  better  information  in  other  areas.  A  smaller  At  and  a 
longer  data  record  provide  more  information  about  the  wheelset  high  frequency 
response  mode  and  low  frequency  model  respectively.  This  increase  in  the 
quantity  of  information  could  be  considered  better  from  the  maximum  likelihood 
processor  point  of  view. 
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Appendix  B 


3P(t. 


SOLUTION  TO 


The  matrix  Riccati  equation  is 


PU)  =  F(t)P(t)  +  P(t)FT(t)  +  G(t)Q(t)G7(t)  -  P(t)H1(t)R~1(t)H(t)P(t) 
Taking  the  partial  derivative  of  both  sides  with  respect  to  0: 


(B.l ) 


3P(t)  _  ,  v  3P(t)  3F(t) 

30  "  Fltj  30_  39_  P 

+  3P£i  F?(t) 

+  G(t)Q(t)  +  G(t)  GT(t)  +  Q(t)GT(t) 

-  p(t)Kv!t)p":(t)H(t)  -  :•  ;i;R?;t)R-1(t) 

-P(t)HT(t)  H(t)P(t)  -  P(t)  ■-1"—  R~1(t^HU)p(t) 

-  KT(t)R'1(t)H(t)P(t)  (3. 

Post-multiply  both  sides  of  (3.2)  by  v.  The  time  notation  will  be  dropped  in 
further  equations  however  all  matrices  are  still  time  dependent. 

f  2.  -  T  §  *  ♦  §  PJL  *  P  Z  *  §  mi  *  ^  IF  Z  *  0  If  ^Z 


*  H  Hr  z  -  PH7"'1  Is  Pi  -  pp" 


3R_1  HPv 


P  1"  R'V  -  f  hTr_1hP^ 


Apply  the  transformations: 


v  _  3P 
*  =  36  2- 


T9 


X  =  39  iL  +  JL 


(B.5) 


3P 


By  reviewing  the  order  of  differentiation  for  the  term  jr.  in  equation(B. 3) 

the  equation  (B.5)  can  be  substituted  in  equation  (B.3). 

X-f  X  =  FX  +  |f 

+  p  w  z  +  If  +  C'Q  W  2-  +  G  ft  gTx  +  -  mV1"! 


~  PhTr_1  f  3S-i  Hib,  -  P  ^  R'VjL  -  §  hV1^ 


(B.6) 


Apply  the  transformation 


•  T  T  -1 
i=-Fjr  +  HR  HP£ 


CB.7)(3) 
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x=Fx+3epx  +  Fifl-x+G^'3try  +  ^^'iy 


ar 

36 


9G 


J  .  3G  „„T 

+  G  -rr  0  y  +  ~ 


36_  -  ’  ”  30  ‘  30  ^  2- 

-1 


-  PHTR-1HX  -  PHTR""1  H  FX  -  pHT  Q  -  P  ~~  R-PHPX 


(B.3) 


Equation  (B.7)  and  (3.8)  are  a  linear  system  and  are  written  in  matrix  form  in 

B.9. 


where 


“ll  W12 


W21  w22 


X 


( £ .  P ) 


KM 


Equation  (B.9)  can  be  solved  using  a  state  transition  matrix  approach. 


y  (t  +At ) 

0 

-p 

< 

o 

-p 

^yY(to’At) 

xCt  +At) 

—  o  _ 

<J)  ( t  ,  At ) 

LyYy  o’ 

V(to’At)- 

.  i(t0L 

where 

$(to,At)  =  exp 

From  equation  (B.lU) 

£(tQ+At)  =  <t)yy(to jAt )  X.(tQ )  +  6yy(to,At)  x(tQ) 

X(tQ+At)  =  <{>Yy(to,At)  jr(tQ)  +  cf>YY(to,At)  ±(tQ) 

Applying  equation  (B.M  and  combining  (l.l6)  and  (1.17): 

9P(t  +At)  9P(t  )  -1 

dB_  =  ^Yy^o’At)  +  99^“  ^ 

9P(t  )  -1 

[6  (to,At)  +  *  (t 0,At)  ~ ge2-] 
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Appendix  C 

DEFINITIONS  OF  VARIABLES  USED  IN  F  AND  G  MATRICES 


Parameter 

Symbol 

Value 

Units 

Wheel set  mass 

M 

12.0U5 

kg 

Nominal  longitudinal  creep  coefficient 

fll 

8705.5 

n 

Nominal  lateral  creep  coefficient 

f22 

25^9.0 

n 

Lateral  spring  constant 

K 

y 

11490.0 

n/r. 

Yaw  spring  constant 

V 

> 

77.3 

n-m,/  rad 

One-half  distance  between  contact  points 

£ 

0.11*7 

n 

Wheel  conicity 

ot 

0.05 

rad 

Wlieelset  rolling  radius 

r 

0 

0.083 

m 

Lateral  damning  constant 

c 

96.5 

n-m-s ;ra 

Forward  velocity 

V 

3.021 

m/s 

Appendix  D 
3P(t.) 

SIMPLIFICATION  OF  — —  FOR  WHEELSET  PROBLEM 

ov 


Simplifying  equations  (B.ll)  through  (B.lM  from  Appendix  5, 
(B.10)  reduces  to: 


T  T  -1 
-F  +  H  R  HP 


,Ep  +  pi: 

L  36  '  S0 


T  -l 

F  -  PH  R  H 


X 


Equation  (D.l)  is  solved  using  the  state  transition  matrix  annroe 


~i_{t  +At  r 

r^(At) 

<5  ( At )  “ 

fv(t  n 

0 

yy 

yY 

0 

Y(t.  +At) 

w  -  O  “ 

-c>Yy(At) 

! 

-p 

<1 

>- 

y 

o 

r 

H 

O 

1  ^ 

$(At)  =  exp 


T  T  -1 
■  -F  +  H  R  HP 


3F.  „  9F" 

Lae  r  96 


F  -  PH1R~^H 


Using  the  same  methodology  as  in  Appendix 


9P(t  +At ) 

_ o _ 

90 


’  -  > 

V{At)  + 


3P(t  )  -1 

[V(At)  +  VAt) 
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Equation  (D.M  is  an  iterative  equation  and  is  solved  for  7-p- 

Of  . 


Because  the  wheelset  is  a  time  invariant  system 
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is  solved 
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for  only  once  for  each  iteration  of  the  maximum 
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Appendix  E 

LISTING  OF  THE  MAXIMUM  LIKELIHOOD  PROGRAM 


no  non 
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CiuL  INDUMP 

THIS  PRO SR AH  IS  THE  CONTROLLING  PROGRAM  FOR  A  SERIES  OF 
SUBROUTINES  1HICH  FOFH  THE  MAXIMUM  1IKELIHOOD  PARAMETER 
IDENTIFICATION  PROCESSOR  FOR  THE  SIMULATED  NHEELSET  DATA. 

COMMON  P  (3 , 3)  ,G  (3 , 1)  , H  (3,  3)  ,  P  (3, 3)  ,  C  ( 1 ,  D  ,  K1F,  K1G,K2G,K1B,  K2H,  K1Q, 
$K 1 d,  RKGM ( 3, 3)  , F1PHIX, K2PHIX ,K 1PHIZ , K2PEIZ, 

*K1aKGM,K2  RKGM,RHASS,P11,P22,FKY,RKSi,RLE,  ALPHA,  RZERO , DT,  NR  ECF. , 
$NRECD,IROSTA, I F OEND, IR 1 STA , IR 1  END, I EATSA ,1 DATEN, NDATA, NRCAL, 

SATTEN  (4,3) ,RVCAI,TVEL,GRCAL,NFOINT,SDl,P (3,3)  ,  B  (3 , 3)  ,  K  IB,  K  IP, 
SRLIKE, K1DPH,K2EPH,K1DRKG,K2DRKG,NPT1,X2C (3,1)  ,BZZRO,  K1DELP, 
SK2DELP, K1 DELK , K2DELK , K 1  DP, K2DF, I FL AG 
DIMENSION  DrLI AP  (3,3) 

DIMENSION  DELIAN  (3,3) 

DIMENSION  DPHI  X  (3,3) 

DIMENSION  DFMTX  (3,3) 

DIMENSION  ZMEAS  (3,1001) 

DIMENSION  ZMEAS1  (3,1001) 

DIMENSION  XHAI  (3,1001) 

DIMENSION  RNU  (3,1001) 

DIMENSION  ZHAT  (3,1001) 

DIMENSION  DXH AT  (3, 1001) 

DIMENSION  XMEAS1  (3,100  1) 

DIMENSION  XMEAS  (3, 1001) 

DIMENSION  XERRCR  (3, 1001) 

DIMENSION  XVECll(ICOI) 

DIMENSION  YVECT 1 (1001) 

DIMENSION  Y VECT2  (1001) 

DIMENSION  U(1,1C01) 

READ (5, 10)  K1F,K1G,K2G, K1H, F2H,K1Q,K  IP 
10  FORMAT(I3) 

WRITE  (6,15)  K1F,R1G,K23,K1H,F2H,K1C,K1F 
15  FORMAT (//, 1  • , 5X, 'K1F= ' ,13, 5X,» K1G='  ,13 , 5X, • K2 G= '  , 13 ,5 X , • K  1H= • , 
$13, 5 X,'  K2H=',13,5X,'K13=*,I3,5X,  *K1E>*#13) 

K 1d*LK=K1 F 
K2JELK=K1 ? 

K1DELP-K1F 
K2JELP=K1F 
K 1 DPH=K1 F 
K2jPH  =  K  IF 
K1  J?=K1F 
K2dF=K1F 

READ (5, 25) ID ATS A, ID ATEN 
25  FORMAT (14) 

NPjINT=  (IDATEN-IDATSA) ♦  1 

NPU=NPOI  NT*  1 

IF±.AG=0 

WRITE (6,45)  IDA1SA,IDATEN,NP0INT,NPT1 
45  FORMAT (//, *  • ,5X,»IDATSA=*,I4,5X, 'IDATiN=  •  ,14,  5X, »NPOINT=', 

$14,  5X, *  NPT1®  * ,1 4) 

CAi.L  DRIVER  (ZME  AS,  ZHE A  S  1,XH  AT ,RNU,ZH  AT,  DELTAP,  DELTAK  , DXHAT, 
$DPdlX,DFMTX,XERBOF,XMEAS1,XMIAS,XVEC11 , Y VECT 1 , Y VECT2 ,0) 

STJ? 

END 


SUBROUTINE  DRIVER (ZHEAS, ZME AS1 ,XHAT,RN U, ZHAT, DELTAP , DELTAK, DXHAT , 
$D PMI X ,D FM  TX , X ERROR ,XMEAS1,XHEAS,XV  ECT1 , YVECT 1 , Y  VECT2, U) 

COMMON  F  (3,3)  ,G  (3,1)  ,H  (3,3)  ,B(3,3)  ,Q(1,1),K1F,K1G,K2G,R1H,K2H,K1Q, 
$K 1 R,RKGH (3,3) , R 1PHIX, R2PHIX , K 1PHIZ , K2PEIZ, 

SKIitKGH,  K2RKGM,  ENASS,  F1 1  ,F22,  FKY ,  RFSI  ,R  IE,  ALPHA,  RZERO  ,  DT,  NRECR, 


'jjfmcinulfi  jfcac-c* 
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SNRaCD, IROST  A,I ROEND, IR 1ST A, I FIEND, IBITSA,IDAT£N, K DAT A, NRC  AL, 
SATXEN  (4,3),RVCAI,TVEL,GBCAL,NP0INT,SET,P(3,3) ,B (3 ,3) , K1B , K IP, 

$R LIKE, K ID PH, K2DFB ,K1 DRUG, K2DRKG, NPT1,  XIC  (3, 1)  , BZEPO, K1DELP, 
SK2JELP, K1DELK, R2DELK, K1£F,K2DF,IFIAG 
DOJBLE  FRECISIC  N  PBIZ  (6,6) 

DIMENSION  PHI X  (3,3) 

DOJBLE  PRECISION  PBIXDP(3,3) 

DIMENSION  ZNEA  S  (3, NPOI NT) 

DIMENSION  ZMEAS1  (3,NPT  1) 

DIMENSION  DXH  AT  (3 , NPOI NT) 

DIMENSION  RNU  (3,NPOINT) 

DIMENSION  ZBAT  (3  ,  NPOINT) 

DIMENSION  DELTAP(K1DELP,K2DELP) 

DIMENSION  DELTA* (X 1DELK , K2DELK) 

DIMENSION  DPHIX  (K  1 DPH  ,  K2DPH) 

DIMENSION  DFMTX  (K1DF,K2DF) 

DIMENSION  Z (6,6) 

DIMENSION  XHAT  (3, NPOINT) 

DIMENSION  XMEAS1  (3, KPT  1 ) 

DIMENSION  XHEAS  (3, NPOINT) 

DIMENSION  XERBCR  (3, NPOINT) 

DIMENSION  XVECT1  (NPOINT) 

DIMENSION  YVECT1 (NPOINT) 

DIMENSION  YVECT2 (NPOINI) 

DOJBLE  PRECISIC N  PH1K(3,3) 

DIMENSION  I XI  T  1  ( 1 0) 

DIMENSION  ITIT2  (10) 

DIMENSION  JLAB1  (10) 

DIMENSION  JLAE2  (10) 

DIMENSION  JLAB3  (10) 

DIMENSION  ITIT3  (10) 

DIMENSION  VD2LJ  (20) 

DIMENSION  VDZLJ2  (20) 

DIMENSION  VRLIKE  (20) 

DIMENSION  V THETA  (2  C) 

DIMENSION  ITIT 5  (10) 

DIMENSION  ITIT4  (10) 

DIMENSION  JLAB4  (10) 

DIMENSION  JLAE5  (10) 

DIMENSION  ITIT6  (10) 

DIMENSION  JLAE6  (10) 

DIMENSION  JLAB7  (10) 

DIMENSION  IIIT7  (10) 

DIMENSION  VCCNST  (2  0) 

DIMENSION  VDETB  (20) 

DIMENSION  U (1 , KPT  1 ) 

TV £L=  3. 02  1 
CALL  PARAH  (TBETA) 

W  RITE  (6, 2  0)  5HASS, F11,?22, SKY, RK  SI  ,RLE  ,  ALPHA  ,R  ZERO,  THETA  ,B  ZERO, 
SDT, SDT 

20  POxtMAT  (// ,  *  •  ,2X,  »RHASS=»  ,F12.5,2X,  •  F 1 1=  »,F12.  5,  2X,  •  F22*  •  ,  F12. 5, 
$2X,*  RRY='  , F 12. 5,2X,'RKSIs',F12.5,2X,'RLEs* ,F12.5, //, 

$•  •  ,  2  X  , *ALPHA= ' ,F1 2.5,2 X, 'RZEFO=' ,F12.5,2X,' THETA«* , 
iFl2.5,2X,»BZ£RC*»,F12.5,2X, 'DT=',F12.5,2X, *SDT= • ,F 12. 5) 

WRITE  (6,  40) 

4  C  FOhMAT(/,'  •,  12X,'G  BATFIX') 

CAi*L  OUTPUT (G, K  1G, K2G) 

WRITE  (6,60) 

60  FORMAT (/, '  * , 12X, 'H  BAIEIX') 

C XuL  OUTPUT  (B, FIB, K2H) 
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WRITE  (6,80) 

80  FORMAT (/, •  *,12X,'Q  MATFIX') 

CALL  OUTPUT  (Q,F1Q,K1Q) 

WRj.TR  (6,100) 

100  FORMAT (/, '  ', 12X,'R  MATRIX' ) 

CAi.L  OUTPUT  (R, KIR, KIR) 

DO  120  I® 1, NPT 1 

READ  (5,  140)  (ZMEAS1  (J,I)  ,J=1  ,3) 

120  CONTINUE 
140  FORMAT  (3P  15. 8) 

DO  160  J®1 ,3 
I=IDA  TSA- 1 
XIZ(J,1)=ZMEAS1  (J,I) 

160  CONTINUE 

WRITE  (6,1  80) 

180  FORMAT (/, '  *,12X,*XIC*) 

CALL  OUTPUT  (XIC,3, 1) 

DO  200  J=  1,  3 

DO  200  1= ID ATS A  , ID ATE N 

M=I-1 

ZM£ AS  (J  ,M )  ® ZM EAS 1  (J,I) 

200  COaTlNUE 

READ  (5,  210)11111 
210  FORMAT  (10A4) 

READ  (  5,  210)  JLAE1 
READ(5,21 0) ITIT2 
READ ( 5,  210)  JLAE2 
READ  (5, 21 0) ITIT3 
READ  (  5,  210)  JLA33 
REhD  (5,210)  IT  IT  4 
REA D  ( 5,  21  0)  JLAB4 
READ  (5,210) ITI15 
S EA D ( 5 , 2 1  0)  JLAE5 
REaD  ( 5,  210)  ITIT6 
R£AD(5,21  0)  JLAB6 
READ  ( 5,  21  0)  I1I17 
READ (5, 21 0) JLAE7 
DO  212  1, NPT 1 

212  READ  (5,213)  0(1, J) 

213  FORMAT  (El  5.8) 

DO  800  ISTEP® 1  ,20 
H RaTE  (6,215)  ISTEP 

215  FORM AT (////, *  •  , 1 X, ' ISTEP®' ,13) 

WRITE  (6,2 1 6) THETA 

216  FORM  AT (/, '  ',  2X,'THETA=*,F15.8) 

CAi.L  SYSNTX (THETA) 

WRITE  (6,2  20) 

220  FORMAT (/, *  ',12X,'F  MA1FIX') 

CALL  OUTPUT (F,K1F,K1F) 

CALL  ZMTBX (Z, K 1Z) 

WRITE (6, 240) 

240  FORM  AT (/, *  »,12X,'Z  MATRIX') 

CALL  0UTPUT(Z,K1Z,K1Z) 

CALL  STMZ  (Z,DT,PHIZ,K1Z,K 1PHIZ , K2PH1 Z) 

IF  (I FLAG  .EQ.  1)G0  TO  1000 
WRITE  (6,260) 

260  FORMAT  (/,  '  12X, 'KALMAN  FILTER  STM  (PHIZ)') 

CALL  OUTDP (PHIZ, 6, 6) 

CALL  R1CATI (PHIZ, K 1PHIZ , P, K 1 P) 

IF(IFLAG  . EQ.  1)GO  TO  1000 
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W  RITE  (6,280) 

280  FORMAT  (/,  1  ',12X, 'STATE  COVARIANCE  RAT  RI X  (P)  •) 

CA  i.L  OUTPOT  (E,  K1P,K1P) 

CALL  KGHTRX 
WRITE  (6,300) 

300  FORMAT (/, '  ' ,121, 'KALMAN  GAIN  MATRIX  (FKGM) ' ) 

CALL  OUTPOT (RKGB,K1FKGM ,K2RKGM) 

CALL  STHF  (F, DT  ,PHI  XDP ,  R  1i  ,  K 1  PHI X  ,K  2PHI X) 

CAi.L  DPTOSP  (PHIXDP,K1PHIX,K2PHIX,PHIX) 

WRITE  (6,320) 

320  FORMAT(/,'  ',  12X,  •  SYSTEM  STATE  TRANSITION  MATRIX  (PHIX)  • ) 

CA-L  OUTPOT  (PHIX,  R  1PHIX ,R2PHIX) 

CALL  STATE  (ZMEA S,XHAT, PHIK,  K1PHIK,  U) 

339  CALL  RINOVA  (XH AT ,ZMEAS , ZHAT , R NU) 

CALL  ZCOV 

WRITE  (6,340) 

340  FORMAT  (/, *  • , 1 2X ,  • MEASU REMENT  COVARIANCE  MATRIX  (B)') 

CALL  OUTPUT  (B,3,3) 

CA-L  RMLPI ( RNU , SC, SD) 

VCJNS2 (ISTEP) =SC 
VDET3  (ISTEP)  =  SD 
WRITE (6, 360) 

360  PORMAT (/, '  ',12X, 'VALUE  CF  THE  LIKELIHOOD  FUNCTION  (RL1KE) •) 

WRITE  (6,380)  R1IKE 

380  FORMAT  ('  ' , 12 X , ' RLIK E= • , S 15 . 8) 

CA*.L  P A RF (THETA, DFMTX) 

WRITE  (6,400) 

400  FORMAT (/,  *  • ,  12X, • FARTIAL  OP  P  IRT  THETA  (DFMTX)  ' ) 

CA-L  OUTPUT (DFPTX,K1DF, K2DP) 

Ckuh  PDP( DFMTX, EE LTAP, CELT A  K) 

WRITE  (6,440) 

44C  FORMAT (/, '  ', 12X, » PARTI AL  OF  P  WRT  THETA  (DELTAP)') 

CAuL  OUTPUT (D2LIAF,K1D31P,K2DELP) 

WRIIE  (6,460) 

460  FORMAT  (/,  •  ',  12X, '  PARTIAL  OF  RKGM  WRT  THETA  (DELTAK)  ') 

CA-L  OUTPUT (DELTAK,K1DELK , K2DELK) 

CAx.L  XSEN  (P3IK,K1PHIK,  DFMTX,  DELTAK,  XH  AT,  ZNEAS,  DXH  AT) 

CAuL  GRADE (EXH AT, RNU,D£LIAP,DELTAJ) 

WRITE  (6,480) 

480  FORMAT (/,'  ' , 12X, ' GPADIENT  -  PARTIAL  CF  THE  LIKELIHOOD  FUNCTION', 
S' WRT  THETA') 

WRITE (6, 500)  EEITAJ 

500  FORMAT  (•  '  ,  12 X  ,  '  DELTA J=  •  , E  1 5.  8) 

CALL  FISHER (P NU , DXH AT, E EL TAP, EELJ2) 

WRITE  (6,520) 

520  FORMAT(/,  •  ',  12X, 'FISHER  INFORMATION  MATRIX') 

WRITE  (6,540) DEI 02 

540  FORMAT  ('  ',12X,'DELJ2=',E15.8) 

CALL  STEP (T  HEX A .DELTA J , DE LJ 2 , DTH ET A , TH ET AN) 

WRITE(6,S50)  DTHITA 

550  FORMAT(//,»  • ,12X,'DTHBTA=* ,E15.8) 

WRITE ( 6, 560) IHETAN 

560  FORMAT(/,'  »  ,  1 2X, ' THET A N= • , E 1 5. 8) 

VTriETA (ISTEP) *THETA 
VDELJ  (ISTEP)  =  D£LTA J 
VDRLJ2  (ISTEP)  *EILJ2 
VR*.IKE  (ISTEP)  »RIIKI 
THi-TA*THETA*0.05 

800  CONTINUE 

805  CONTINUE 


on  o  n  o  o 
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HRITE  (6,8 10) ITIT3 
810  FORMAT ('1',40X,10A4) 

CALL  HPL0T1 (VTHITA, VBLIKE,ISTEP#4  0, JIAE3) 
HRITE  (6,81 0) I1IT1 

CALL  HPLOT1  (VIHIT  A  ,  VCON  £X,  I  STEP,  4C,  JLAE1 ) 
HRiTE  (6,8  10)  IT1T6 

CALL  HPLOT1 (V  THETA, VEETB, ISTfcP, 40 , JLAE6) 
MRITE  (6,8  10 )ITIT4 

CAi,L  MPLOT1  (V1HETA,VDELJ,ISIEP,40, JLAB4) 
HRITE  (6  ,8  10)  IT1T5 

CALL  HPLOT1 (VTH ETA , VDELJ2 , 1  STEP ,4 0, Ji A ES ) 
900  CONTINUE 
1C0C  RETURN 
ENJ 


SU3POUTIN E  PAEAM (1UETA) 

SUBROUTINE  PARAM  READS  IN  THE  PARAMETERS  NECESSARY  TO  FORM 
THE  SYSTEM  P A1EIX 

COMMON  F(3,3)  ,G  (3,  1)  ,H  (3,3)  ,R  (3,3)  ,Q  (1,1)  ,!C  1 F, K 10,  K2G,  K1H,  K2H,K1Q, 
SK1r,RK3M  (3,3)  , K1 PH1X , K 2FHII , K 1  PHIZ, K2PKIZ, 

$K1RKGM,K2RKGM ,RMASS,F1 1  , P22  , R  KY , R  FSI , RLE,  ALPHA,  RZ  ERO,  DT,  KPECF., 
IN3£.CD,IR0STA,IB0END,IR1STA,  IR1END,1DATSA,IDATEN,  NDATA,NPCAL, 

SAT  TEN  (4,3)  ,  R  VC  A I  ,T  VEL  ,  GRC  AL  ,N  POINT  ,S  ET  ,P  (3  ,3)  ,B  (3 , 3)  ,  K1E,K  IP, 
$FLIKE,K1DPH,K2DFH,K1DRKG,K2DRKG,NPT1,XIC (3,1), BZERO, K1DELP, 

SK  2jELP, K1DELK, K2DELK , K 1  DP , K 2 DF, I F LAG 
1C  FORMAT  (FI  5.  8) 

PE AD (5,  10)  5HASS,F11,F22,RKY,RKSI,RLE,  ALPHA,  PZ  ERO,  THETA,  BZERO,  DT, 
$S  DT 

DO  20  1= 1 , K  1G 
DO  20  J=1,K2G 
20  R  E^D (  5,  10)  G  ( I ,  J) 

DO  30  1=1, K1H 
DO  30  J=1,K2H 
30  READ  (5,  10  )  d(I,J) 

DO  40  1=1, K1Q 
DO  40  J=1,K1Q 
40  READ  (5,  10)  Q(I,J) 

DO  50  1=1, KIR 
DO  50  J=1 , K 1 R 
50  READ(5,  10)  P(I,J) 

DO  70  1=1,4 
DO  70  J=1, 3 
REhD  (5  ,  10 )  ATTEN  (I  ,J ) 

70  CONTINUE 
RETURN 
EN  J 


SUBROUTINE  SI SFTX  (THETA) 

C  SUBROUTINE  SISMTX  FORMS  THE  SYSTEM  MATRIX  -  F  MATRIX 

COMMON  F  (3,  3)  , G  (3,1)  , H  (3, 3)  , R  (3,  3)  ,  Q  ( 1 ,  1)  ,K1F,K13,K2G,K1B,  K2H,  K  la, 
SK  1  d,  RKGH  ( 3,  3)  ,  R1PHIX,K2PHIX,K  1  PHIZ  , K2PHIZ, 

SK 1dKGM,E2RKGM,RPASS, F1 1 ,F22, RKY , RK SI , RLE , ALPHA, RZ ERO , DT, NR  BCR, 
SNRiCD, IRO STA, IR  CEN D, IR 1ST A, IR 1  END, IDATSA , IDA TEN, N  DATA, NRCAL, 

SATTEN  (4,3) , RVCAL,TVEL,GFCAL, NPCI NT, SDT, P (3,3) , B (3,3) ,K1B,K IP, 
SRLIKE,K1DPH,K2EPH,K1DRKG, K2DEKG , NPT 1 ,X IC  (3,1) , BZERO, K1DELP, 
SK2dELP, K1 DELK , K2DELK, K 1  DP, K  2D F, I FLAG 
F ( 1 ,  1)=-(  ((  (2.C)*F22*THETA)  *  ( BZERO *TV El)  )  /(RHASS*TVEL)  ) 

F  (1,2)  =2.0*F22MHETA/RMASS 
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P  0.3)*  (-aKD/RHASS 
F  U,1)*0.0 

F(2,2)  =  (-RKSI)*1VEL/(2.0*F11*THETA*  (RLI**2.0)  ) 

P  (2 ,3)  =  (•  ALPHA )  *TVEL/  ( B  IE  *BZER0) 

F(3,1)  =  1.0 
F  (3,2)  =0.  0 
F(3,3)=0.0 
REIUFN 

ENu  # 

C 

C 

SUBROUTINE  ZMTFX(Z,K1Z) 

C  SUBROUTINE  2KTRX  FORMS  THE  Z  MATRIX  FOR  THE  KALMAN  ENGLER 

C  SOLUTION  TO  THE  R1CCATI  EQUATION 

COMKON  F(3,3)  ,G  ( 3,  1)  ,H  (3, 3)  ,  R  (3  ,3 )  ,0  (1  , 1  )  , K1  F , K1 G  ,K2G,  K1 H,  K2H  ,  K IQ  , 
SK In, R KGM  (3,3) , F 1 PHIX , F 2 FHIX , K tPHIZ , K2PHI Z, 

EKlBKGM,  K2RKGM,  EMASS,F1  1  ,P22  ,P.FY,RKSI  ,R  IE, ALPHA  , F ZERO  ,DT,  NFECF  , 
iN  P  fiCD , IPO ST A, I POEND, I F 1 STA , I R 1 END ,1 C  AT S  A, I D ATEN ,NCATA,NFCAL, 

$A TIEN  (4,3)  ,  RVCAL,TVEL,GFCAL  ,N POINT,  SET  ,P  (3,3)  ,B  (3 ,3)  ,K16  ,K1P, 
*RLIKE,K1DPH,K2CPH,K1DFKG,K2DEKG,NFT1,X1C  (3,1)  ,BZERO, K 1 DLLP , 
$K2DnLP, K1DELK,K2DELK,K 1 EP , K2 DF, I  FLAG 
DIMENSION  Z  (6 , 6) 

DIMENSION  FI (3, 3) 

DIMENSION  Z11  (3,3) 

DIMENSION  FI  (3,3) 

DIMENSION  HI ( 3 , 3) 

DIMENSION  FI H  (3,3) 

DIMENSION  1 12  (3,3) 

DIMENSION  GT  (1,3) 

DIMENSION  OGT  (3,3) 

DIMENSION  Z2  1  (3,3) 

DIMENSION  Z22  (3  ,3) 

DOUBLE  PRECISION  VKAREA  (9) 

DOu  3L  E  PP.ECISIC  N  FIDP(3,3) 

DOUBLE  PRECISION  PDP(3,3) 

01  CA*.L  TRANS  (F,  FT ,  F  IF ,  K  1  f  ,  K 1  FT  , K2FT ) 

CA LL  N2G(FT,Z11,K1FT,K2FT,K1Z11,K2Z1  1) 

IDGT=0 

CALL  SPTODP  (F,K1R,K1S,REP) 

CA..L  LIN  V  IF  (KDP,K1R,K1R,RIDP,IDGT,NKARIA,IFP) 

CALL  DPTOSP (RIEF,K1F,K  IF, FI) 

CA«L  TRANS(H, HT , K 1 H, K2 H , K 1HT , K2HT) 

CAi*L  MULT  (RI,H,RIH,K1R,K1F,  K  1H,  K2H,  K  1FIH,  K2RIH) 

CAwL  MULT (HI,F1H,Z12,K 1  HI , K 2HT, K 1  FI H ,K2RIH, K1 Z 12 , K2Z 12) 

C A^L  TRANS(G,GT,F1G,K25,K1GT,K2GT) 

CALL  NULT(Q,GT,QGT,K1Q,K1Q,K1GT,K2GT ,K  1QGT,K2 QGT) 

CAnL  MULT  (G,QGT,Z21,K1G,K2G,K1QGT,K2CGT,K1Z21,K2Z21) 

CA^L  SCALAR  (F  ,  1 . 0 ,Z22  ,  K  IF  , K  IF,  K 1Z22  ,  F2Z22) 

DO  10  1=1 , K 1 Z1 1 
DO  10  J=1,K2Z11 
B*1 
N=J 

10  Z  (S,N)=Z11  (I,J) 

DO  20  1=1 , KlZ 12 
DO  20  J=1,K2Z1 2 
M  =  I 

N=j  *K  2Z  1 1 

20  2  (M,N) =Z12  (I, J) 

DO  30  I*1,KlZ2l 
DO  30  J=1,R2Z21 
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M*I*K1Z11 

N=J 

Z  (d# N )  - Z2  1  II,  J) 
DO  40  1=1 ,K1Z22 
DO  40  J=1,K2Z22 
K=i+K1Z21 
N=J»K2Z12 
Z  (M,N)=Z22(I,  J) 
K1Z=K 1Z 11+K1Z2  1 
RETURN 
END 


SUBROUTINE  STMf  <F,CT,PH1,K1F,K1PHI,K2PHI) 

SU  BROUTIN2  STM  F  DIMENSIONS  THE  NEC ESS  EFT  ARRAYS  FOR 
SUBROUTINE  STM.  IT  ACTS  AS  A  DRIVER  FOR  SUBROUTINE  STM. 
DOUBLE  PRECISION  FDP(3,3) 

DIMENSION  F (K 1 F,K 1 F) 

DOUBLE  PRECISION  PHI(K1F,K1F) 

DOUBLE  PPBCISICN  FDT  (3, 3) 

DOUBLE  PRECISION  PRD  (3,3) 

DOUBLE  PRECISION  PRD1(3,3) 


D0U3LE  PPBCISICN  FDT  (3, 3) 

DOUBLE  PRECISION  PRD  (3,3) 

DOUBLE  PRECISION  PRD1(3,3) 

DOUBLE  PRECISION  TERM  (3,3) 

DOUBLE  PRECISICN  FHI1(3,3) 

DOUBLE  PRECISICN  EIFF(3,3) 

DOUBLE  PRECISION  RM(3,3) 

C  A  uL  SPT0DP|P,K1F,K1F,FD?) 

CA^L  STM  ( PDP, PHI, FDT,PRD,P2D1, TERM, PHI  1,  DIFP,  RM,  DT,  K  IF, 
iKlPHI,K2PHI,IFLAG) 

RETURN 

tN  j 


SUBROUTINE  STMZ  (Z , DT , PHI, K 1Z , K 1  PHIZ , K2 EH IZ) 

SUBROUTINE  S1HZ  DIMENSIONS  THE  NECESSARY  ARRAYS  FOR  SUB¬ 
ROUTINE  STM.  IT  ACTS  AS  A  DPI VEF  FCP  SUBROUTINE  STM 
DIMENSION  Z  (K  1  Z,K  1  Z) 

DOUBLE  PRECISION  PHI  (6,6) 

DOUBLE  PRECISICN  FDT(6,6) 

DOUBLE  PRECISION  PHD  (6,6) 

DOUBLE  PRECISICN  FRD1(6,6) 

DOUBLE  PRECISICN  TERM  (6,6) 

DOUBLt,  PRECISICN  ZDP(6,€) 

DOJ3LE  PRECISICN  PHI1(6,6) 

DOUBLE  PRECISICN  DIPF(6,6) 

DOUBLE  PRECISICN  RM(6,6) 

CAi,L  SPTODP  (Z,R1Z,K1Z,ZD?) 

CALL  STM(ZDP, PHI, FDT, PRC, PRD  1,  TERM,  FBI  1,  DIFF, R M, DT, K 1Z, 
SK1PHIZ,  K2PHIZ,2FLAG) 

RETURN 

ENJ 


SUBROUTINE  STM  (f , PHI , PCI, PHD, PPDl ,TEFM ,P H1 1, DIFF,RH, DT , K IF, 
iK  1?HI , K2PHI , I FLAG) 

SUBROUTINE  STM  CALCULATES  THE  STATE  TRANSITION  MATRIX. 
DOUBLE  PRECISION  F  (K1F,K1F) 

DOUBLE  PRECISICN  FDT  (MF,K1F) 

DOUBLE  PRECISION  PHI(K1F,K1F) 

DOUBLE  PRECISICN  PBD  (R 1  F,  K1F) 
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01 


20 


120 


125 


160 


260 

260 

300 

320 

3**0 

360 

360 

400 

420 

440 

450 

460 


DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOJBLE  PRECISION 
DOUBLE  PRECISION 
DOJBLE  PRECISION 
DTDP=DBLE  (DT) 
CT=DBLE  (1 .01-16) 
IPLAGO=0 


PH1 1  (KIP, KIP) 
PfiDI  ( K  1 P,  K 1  F) 
DIFF  (K  IP,  K  IP) 
TERN  (KIP, KIP) 
Rfl  (KIP, KIP) 

B 

COEF 
CT 
CHK 
DTD  P 
FACTDP 


B=DBLE(0.0) 

CA-L  SCALDP  (F,B,Rfl,KlF,!MP,K1RM,K2RB) 

DO  20  1  =  1 ,  K 1 R  fl 
Rfl  (I,  I)  =D3LE  (1.  C) 

DO  320  1=1,200 

IF  ( I  .GT.  1) GO  TO  120 

CALL  SC ALDP (F , DTDP , PDT , K 1 F, K IF, K 1 FD1, K 2FDT) 

CA*.L  ADDDP(RK,FLT,FH1,K1RM,K2RH,K1PHI,K2PHI) 

B=DBL£  (1.0) 

CAi.L  SC  ALDP  (P,  E, PRD, KIP, KIP  ,K1PFD,K2PSE) 

GO  TO  300 

CALL  HULTDP(F,FFD,PED1,K1F,K1F,K1PRC,K2PRD,K1PRD1 ,K2PRD1) 
CHK=DELE(1.0E  71) 

CKuL  0VERDP(PRC1,CHK,K1FRD1 ,K2PRD1 ,1FL AGO) 

IF  (IPLAGO  .  EQ.  0)GO  TO  160 
WRITE  (6,125) 

FORMAT  (/, '  '  ,2X,'?F.D1  EXCEEDS  EXPONENT  OVERFLOW.', 

$1X, 'CURRENT  VALUE  OF  PHI  USED.  EXECUTION  CONTINUES.') 

GO  TO  U  20 


B  =  J3LL  (1.0) 

CA-L  SCALDP  (PRE1,B,PaD,MPaD1,K2PFE1,R  1PRD,K2?RD) 

COEF=  (DT*  *1) / ( FACTDP (I )  ) 

CHK=D BLE  (1. OE-76) 

IF  (CO EF  .  LE  .  CHK)  GO  TO  380 

CA*.L  SCALDP  (PRE,C0IF,TERH,K1PRD,K2PRE,K1  TERM,  K2TERK) 

CALL  ADDDP (PHI .TERM , PHI ,K 1PHI, K2PHI, K 1PHI, K2PHI) 

CAuL  SUBDP(PHI,PHI1,DIPF,K1PHI, K2PHI,K 1DIFF, K2DIPF) 

DO  260  J=  1 ,  K1  PHI 
DO  260  K=1,K2PHI 

IF  (DABS  (DIFP(J,K)  )  .GT.  CT)GO  TO  280 

CONTINUE 

GO  TO  420 

IP(I-200)  300,340,340 
B=wBLE  (1.0) 

CAkL  SCALDP (PHI, B, PHI 1,K1 PHI ,K2PHI , K 1PHI 1 , K2P HI 1) 

CONTINUE 
WRaTE  (6,360) 

FORMAT (//, '  ' , 5X, ' 200  ITERATIONS.  NC  CONVERGENCE* ) 

GO  TO  450 
V RITE(6,  400) 

FORMAT  (//,'  '  ,2X,  '  COEF  .LE.  1. 0E-76.  ITERATION  STOPPED.', 
SIX, 'CURRENT  VALUE  OP  PHI  USED.  EXECUTION  CONTINUES.') 
WRITE  (6,440)  I 

FORMAT {//, '  ', 2X, 13, IX, 'ITERATIONS  POR  STB.') 

WRITE  (6,460) 

PORRAT (/, '  ' , 12X, ' STB  DIFF  MATRIX') 
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CA-L  OUTDP  (DI  FF,K1 01 FF,  K2DIPF) 

480  BEIURN 
END 
C 
C 

SUBROUTINE  RICATI  (PHIZ  ,  K1  PHI ,  PSP,  K  1 P) 

C  SUBROUTINE  FICCATI  CALCULATES  THE  STATE  COVARIANCE  MATRIX  I.E. 

C  THE  P  MATRIX. 

DOJBLE  PR  EC I SIC N  FHI11(3,3) 

DOUBLE  PRECISION  PHI12(3,3) 

DOUBLE  PRECISION  FHI21(3,3) 

DOJBLE  PRECISION  PHI22(3,3) 

DOUBLE  PRECISION  W  (3,3) 

DOJELE  PRECISION  X(3,3) 

DOUBLE  PRECISION  1(3,3) 

DOJBLE  PRECISION  Z(3,3) 

DOJBLE  PRECISION  DIFF{3,3) 

DOJBLE  PRECISION  E  1  (3 , 3 ) 

DOJBLE  PRECISION  P2(3,3) 

DOJBLE  PRECISION  F2T{3,3) 

DOJBLE  PRECISION  ZI{3#3) 

DOJBLE  PRFCISICN  PHIZ (K 1  PHI , K 1 PHI) 

DOJBLE  PRECISION  PSUH(3,3) 

DOJBLE  PRECISION  F(3,3) 

DIMENSION  PSP  (3,  3) 

DOJBLE  PRECISION  CONST 
DOJBLE  PRECISION  ERROR 
DOJBLE  PRECISION  HKAREA  (36) 

DOJBLE  PRECISION  A,9 
01  ERROE=DELE( 1.0E-16) 

Kl l=K1PHI/2.0 
DO  20  1=1, K 1 1 
DO  20  J=  1 ,  K  1 1 
20  PHxll  (I,J)=PHIZ(I,J) 

DO  40  1=1 ,K11 
DO  40  J= 1 , K 1 1 
M=I 

N=S1 1 ♦ J 

40  PHI12  (I,J)  =  PHIZ  (K,N) 

DO  60  1=1, K11 
DO  60  J= 1 , K  1 1 
K=j.  ♦  K  1 1 
N=J 

60  PHI21  (I,  J)  =  ?HIZ  (M,  N) 

DO  80  1=1,  K 1 1 
DO  80  J=1 ,K 1 1 
H=K11*I 
N=K1 1*J 

80  PHI22  (I  ,  J)  =PHIZ  (M,  N) 

A*D3LE  (0.  0) 

CALL  SCALDP  (PH1 1 1 ,  A, P2 ,  K1 1,  K  11,  K  IP 2,  K2P2) 

ICJUNT- 1 
85  B=DBLE  (1.0) 

CA..L  SCALDP(P2,E,P1,K1?2,K2P2,K1P1,K2P1) 

CALL  BULTDP(PHI22,P2,H,N11,K11,K1P2,K2F2,K1W,K2W) 

C A«L  ADDDP(PHIi1,B,X,K1 1 , Kl 1 , K IX, K2X) 

CALL  MULT  DP  (PHI12,P2,I,K11,K11,K1P2,K2P2,K1Y,K2Y) 

CAi.L  ADDDP  (PHI11,Y,Z,K1  1  ,K1 1,K1Z,K2Z) 

IDOT»0 

CALL  LINV1P(Z#K1ZfKl1, Z1,IDGT,UK ARE A,I IR ) 


CKuL  HOLTDP (X,ZI,F2,K1X,K2X,K1Z,K2Z, *1P2,K2P2) 

CAi,L  TRAN  SD  (P2  ,P2T ,K  IP 2 ,K2P  2 , K1  P2T ,  K2P2T) 

CALL  ADDDP(P2,P2T,PSUH,K1P2,K2P2, K1PSUM, K2PSOM) 

CON  ST=D BLE ( . 500) 

CALL  SCALDP (PSD E, CONST , E2,K 1 PSOH, K 2 PSD B, K1P2.K2P2) 

CAi,L  SUBDP(P2,P1,DIFP,K1P 2,  K2P2, K1DIFF.K2DIF?) 

DO  100  1=1 , R1DIFF 
DO  100  J=1,K1CIFF 

9C  IF  (DABS  (DIF P  (I,J) )  -ERROB)  100.  100,  1  10 
100  CONTINUE 
GO  TO  150 

110  ICOUNT=ICOUNT*  1 

IF ( ICOONT  .  3T.  200)  GO  TO  115 
GO  TO  85 

115  N  RITE  (6,120) 

120  FORMAT (//, '  ',5X,'200  ITERATIONS.  RICCATI  SOLUTION  DOES  NOT', 

5*  JCN  VERGE.*) 

B  =  J3L£  (1.0) 

150  C A_L  SCALDP (P  2,  B ,  P,  K  1  F 2  ,  K2P2  ,  K  1 P , K  2  F ) 

CA  l,L  DPTOS?(P,K  IP, KIP, PSP) 

IF (ICOUNT  .GE.  200)  GO  TO  165 
SPITE  (6, 160) ICCUNT 

160  FORMAT  (/,  '  • , 5X ,13, 1 X, • ITER ATIONS  FOR  FICCATI  SOLUTION  TO', 

J>'  CON  V EP3  E.  ' ) 

165  SRi.TE  (6,170) 

170  FOxsMAT(//,'  ',  12X, 'FICCATI  DIFFERENCE  PATPIX.') 

C A..L  OUTDP (DIFF.K1DIFF, B2DIFF) 

LETUP N 
ENJ 
C 
C 

SUBROUTINE  KGMTRX 

C  SUBROUTINE  KG  HI  RX  CALCULATES  THE  KALMAN  GAIN  MATRIX. 

COMMON  F(3,  3)  ,G(3,  1)  ,H  (3,3)  ,R  (3,3)  ,Q  (1  ,1)  ,K1F,K1G,K2G,  K1H,  K2H,K1Q 
JKU.RKGM  (3,3)  ,  R 1 PHIX  ,  K  2 IHIX  ,  K  1  PHIZ ,  K2PF.IZ, 

SK  IaKGM,  K2RKGM,PMASS,F1 1  ,F 22, BKY,PKSI, RIB, ALPHA,  3ZE30,  DT,  NF.ECR, 
SN?SCD,IR0STA,IR0ENC,IR1STA,IR1END, ID  AT S A , ID ATEN , N DAT A , NR  CAL, 
SATTBN  (4,3)  ,RVCAL,TVEL,  SRC  A  L,  N  POINT  ,SDT,P  (3,3)  ,  B  (3 , 3)  ,  K  IB  ,  K  IP, 
SRLIKE,K1DPd,K2DPH,KlDPKG,K2DEKG,NPTl,XIC  (3,1), BZEPO, K1DELP, 
£K2J.iLP,  KlDEtK,  * 2EF  I K ,  K  1  DP  ,  K  2DF ,  I PLA  G 
DIMENSION  HTF I  (3,3) 

DIMENSION  H T  ( 3, 3) 

DIMENSION  RI  (3,3) 

D0J3LE  PRECISION  WKAREA  (9) 

DOUBLE  PRECISION  RIDP(3,3) 

DO  J  BLi  PPECISIC  N  RDP(3,3) 

01  IDST=0 

CALL  SPTODP  (R,F1R,K1R,F£P) 

CAi.L  LINY  IF  (RDF,K1P,K1R,RIDP,IDGT,MKARIA,IER) 

CALL  DPTOSP  (R1DF,K1R,K1R,RI) 

CA  i.L  TFANS(H,HT,K1H,K2H,K1HT,K2HT) 

CA..L  MULT  (HT,RI,HTRI,KlHT,K2HT,K1FrK1R,K1HTRI,K2HTRI) 

CA  i.L  MULT(P,HTFI,RKGM,K1P,K1P,K1H1RI,K2HTRI,K1RKGH,K2RKGM) 

RETURN 

LNJ 

C 

C 

SUBROUTINE  STAELE 

COMMON  F(3,3)  ,G  (3,1)  ,R  (3,3)  ,R  (3,3)  ,Q  (1  ,  1)  ,K1F,K1G,K2G,  K1H,  K2H.K1Q 
JK 1 2 , R  KGH (3,3) , K 1PH1X ,K2PHI X,K 1PHIZ, K2P BIZ, 


9*4 


SKUKGM,K2RKGH,FEASS,  PI  1,F22,  BKY , RKSI ,RLE ,  ALPH A, RZERO , DT, NR ECR , 
SNRaCD,  IROSTA,  IROEND,  IR 1ST A, IR 1END, IEATSA , IDA  TEN, N DATA, NRCAL, 

* ATT  EN  (4,3)  ,  RVCAL  ,TV  EL,G  RC  AL,  NPOI  NT,  SDT  ,  P  (3,  3)  ,B(3,3)  ,K1B,K1P, 
$RLIKE,K1DPH,K2DPH,K1DRKG,K2DRKG,NPT1,XIC  (3,1) ,BZERO, K1DELP, 
SK2jELP,K1DEIK, R2DELK,K1DF,K2DF,IPLAG 
DIMENSION  R H  (3,3) 

DliENSION  FKMTX  (3  ,3) 

DOJBLE  PRECISION  HKAR  EA  (6) 

COiPLEX*16  EIGVEC  (3,3) 

COMPLEX*  1  6  HDF  (3) 

DOUBLE  PRECISION  PKDP(3,3) 

CALL  BULT(RKGH,H,RH,KlRK3M,K2RKGH  ,R1H,  K2H,K1RH,K2RH) 

CA^L  SUB  (F,RH,FKMTX,K1F .K1F.K1FK, K2FK) 

WPiTE  (6,  10) 

10  FORM  AT (//, '  12X, 'KALMAN  FILTER  SYSTEM  MATRIX  -  (P-KR)  • ) 

CA^L  OUTPOT(FKMTX,K1FK,K2FK) 

IJJB=0 

CA kL  SPT0DP(FKMTX,K1F,K1F,FKDP) 

HPaTE  (6,15)  IER 

15  F0rMAT(//,'  • ,11, 'IER=' ,14) 

WRITE  (6,20) 

20  FOaMAT (//, '  • , 12X, 'EIGENVALUES  OF  (F-KH) * ) 

DO  40  1=1 , K  1 FK 
40  WRITE  (6,6  0)  HD  F(I) 

60  FORMATC  '  ,  14X,2D15.  8) 

DO  80  1= 1 , K  IF K 

IF  (FEAL  (HDP  (I)  )  .IE.  0)  30  TO  120 
80  CONTINUE 

WRITE  (6,100) 

100  FORMAT  (//,»  «,«*****  FILTERED  SYSTEM  IS  UNSTABLE  ****»') 
IF„AG=1 
120  RETURN 
END 


SUBROUTINE  STATE  (ZMEAS , XHAT , PHIK2 , K1 PHIK,U) 

SUBROUTINE  STATE  CALCULATES  THE  STATE  ESTIMATES  OP  THE 
KALMAN  FILTER  ACCORDING  TO  THE  CONTINUOUS  FORM  OF  THE 
KALMAN  FILTER. 

COMMON  F ( 3,  3)  »G  (3,1)  ,H  (3,3)  ,R  (3 , 3)  , Q  (1  , 1 )  ,K1  F ,  K1 G  ,  K2  G ,  K1 H ,  K2H ,  K  1  Q, 
$KlH,FKGM  (3, 3)  ,K1EHIX ,K2PHIX ,K  1PHI Z, K 2FHI Z, 

SK  litKGH,  K2RK  GM,FMASS,F1  1  ,F22  ,  R  F  Y  ,  R  FSI  ,R  IE ,  ALPHA  ,  RZ  ERO,  DT,  NR  EC  P, 
$NRaCD,IR0STA,IF0END,lF1STA,IR1END,IDATSA, ID ATEN, ND AT A, NRCAL, 

SAT  TEN (4,3) , R VC A I ,T VEL , GRC AL , NPOINT , SET , P  (3 , 3 ) , B (3 , 3)  , K  IB , K  IP, 
$PLIKE,K1DPH,K2DPH,K1DPKG, K2DFKG, NPT1,XIC  (3, 1)  ,BZEFO, K1DELP, 
SK2DELP, K1DELK, K2DELK, K 1 DF, K2DF , IP LAG 
DIMENSION  ZMEAS  (3, NPOINT) 

DIMENSION  XHAT  (3,NP0IN1) 

DIMENSION  FB (3 ,3) 

DOJBLE  PRECISION  1KABEA  (36) 

DliENSION  PHIKI  (3,3) 

DIMENSION  HNEG  (3,3) 

DOJBLE  PRECISION  PHIK(3,3) 

DIMENSION  PHIK SP  (3  ,3 ) 

DIMENSION  FKP  (3, 3) 

DIMENSION  FKF1  (3,3) 

DIMENSION  RM  (3,3) 

DliENSION  T2  (3,3) 

DOJBLE  PRECISION  PHIK2(3,3) 

DOJBLE  PRECISION  FACTOR 


mm 


* 
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DIMENSION  PRD  1  (3.3) 

DIMENSION  PP.D2  (3,3) 

DIMENSION  XM1  (3,1) 

DIMENSION  XE 1  (3,1) 

DIMENSION  XE2  (3,1) 

DIMENSION  GAMMA  (3,3) 

DIMENSION  XI  (3,1) 

DIMENSION  X2  (3,  1) 

DOUBLE  PRECISION  FKFDP(3,3) 

DOUBLE  PRECISION  fKFIDP(3,3) 

DOUBLE  PRECISION  PHIKID(3,3) 

DIMENSION  U (1 ,NPT1) 

DIMENSION  01(1,1) 

DIMENSION  A  1  (3,1) 

DIMENSION  A2  (3 , 1 ) 

DIMENSION  LAMEDA  (3,1) 

DIMENSION  X  3  ( 3 ,  1) 

DIMENSION  XEP  (3,1) 

CA^L  NEG(H,HNEG,K1H,K2H,K 1BNEG, K2HN2G) 

CALL  MULT(RKGM,HNEG,RH,K1RKGM,K2RKGM,K1HNZG,  K2HNEG,  K1RH,  K2F.H) 

CA^L  ADD(F,  RH,  FKP.K1F,  K  1P,K  1FKF,  K2FKF) 

WHITE (6,5) 

5  PORK  AT (//, *  ', 12X, 'KALMAN  PILTEP  SYETEE  MATRIX  (F  KF)  ' ) 

CALL  OUTPUT  (FKF,K1FKF,K2PKP) 

CA*.L  STMF(?KF,D1,EHIK,K1FKP,K1PHIK,K2PHIK) 

PA  C 10  R=  1 . OD  OC 

CALL  SCALDP (PHI K , FACTOR , FHI K2 , K 1  PHI K , K 1PHIK, K 1 PHK 2, K 2PHK 2) 

CALL  DPTOSP (PHIK,K1PHIK,K2PH1K,PHIKSP) 

WRITE  (6 ,10) 

10  FORMA  T ( //, '  ' ,  1 2X, 'STM  FOR  KALMAN  PILTER  SYSTEM  (PHIK)') 

CALL  OUTDP  (PHIK,K1PHIK,K2FHIK) 

I DJT=0 

CAi.L  LINV1F  (PHIK,K1PHIX,K1?HI  K , PHI KI D ,1DGT, W K  AREA,  IE R) 

CA  DPTOSP  (PHIKIE,KlPdIK,K2PHIK,  PHIKI) 

ID  GT=0 

CP.uL  SPTODP  (FKF,K1FKF,K2PKF,FKPDP) 

CALL  LIN  V  IF  (F  KF  DP  ,  KI  FK  F  ,K  1FKF  ,  FKFICP  ,IDGT,  K  K  AP.  E  A,  I  EE) 

CA  UL  DPTOSP (PKFIDP,K1FKF,K2FKF, FKFI) 

CAlL  SCALAR (F,0.0,RS,K1F,K1P,K1RM,K2RM) 

DO  20  1=1 ,K1FK 
FMlI,I)=1.0 
2  C  CONTINUE 

CA i.L  SUB  (RK,PHIKI,T2,K1RM,K2RK,K1T2,K212) 

CALL  MULT  (T2  ,  R  KG  M  ,  PED  1  ,  K  1T2  ,  K  2T2,  K  1  RKG  E,  K2EKGM  ,  K  1  PRD  1 ,  K  2Pn  D 1 ) 

C A  i,L  MULT  (FKPI ,  FRC1 ,  PRD2.K1  FKF,K2FKF,K  1PRD1,  K2PRD1 ,  K  1PRD2,  K2PRD2) 
C AmL  MULT (PHIKSF,PPD2,GAEMA , K 1 PHI K , K 2PHI K, K 1 PRD2 , K2PRD2, K 1G AM , 
$K2GA M) 

CA«,L  MULT  (T2,G,A1,K1T2,K2T2,K1G,K2G,K1A1,K2A1) 

CALL  MULT (FKFI ,A1, A2, K 1 FKF, K2FKF, K 1A  1, F2A1.K1 A2.K2A2) 

CALL  MULT  (PHI  KSP,  A  2,  LA BOA  ,  K  1  PHIK ,  K  1PHIK  ,  K  1A2,K2A2,K1LANB,  K2LAMB) 
CALL  SCALAR(XIC,1.0,XE1,3,1,K1XE1,K2XE1) 

DO  1003  1= 1 .NPOINT 

DO  100  J=1 ,3 

XM  1  (J,  1 )  =  ZM EAS  (J,I) 

100  CONTINUE 

U1(1,  1)  =u  (1,1) 

CA*L  HOLT  (PHI  KSP, XE1,X1,K1  PHIK,  K2PHIK,  K  1  XE  1,  K2XE  1, K 1 X  1.K2X  1) 

CAlL  MULT(GAMMA,XM1,X2,K1GAM,K2GAM,3,1 ,K1X2, K2X2) 

CALL  MULT (L AHEDA ,0 1 ,X3 , K 1LA MB, K2L A  SB, 1 , 1.K1X3.K2X3) 

C A»L  ADD(X1,X2,XEP,K1X1,K2X1,K1XEP,K2XIP) 
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CkuL  ADD(XEP, X3,XE2,K1XEP,R2XEP,K1XE2, F2XE2) 

DO  200  J=  1, 3 
XHAT  (J,I)  =X£2  (J,1) 

200  C03TIN0E 

CAU  SCALAF  (XE2 , 1 . 0,X E 1  ,K  IX  E2 ,  K2XE2 ,  K1XE1,K2X2  1) 

1000  CO ii TIN 0 E 
RETURN 
END 
C 

c 

SUBROUTINE  ZCOV 

C  SOBROUTINE  ZCOV  CALCOALTES  THE  P  MATRIX  FOR  THE  MAXIMUM 

C  LIKELIHOOD  PARAMETER  IDENTIFICATION  PROCESSOR. 

COMMON  F(3,3)  ,G  (3,1)  ,H  (3,3)  ,R  (3,3)  ,Q  (1,1)  , K 1 F,  K 1G,  K2G,  K1 H,  K2H,  K1  Q, 
*K1i,FKGK  (3,3)  , K1 PHIX, K2FHIX , K 1PHIZ, K2PHIZ, 

$K1rKGM,K2RKG3,RMASS,F1 1,F22,RKY ,RKSI  ,R IE, ALPHA, RZEBO, DT, NRECR, 
*NRr;CD,IF0STA,IF0END,IRlSTA,IR1END, ID  AT  £ A, ID ATEN, ND AT  A, NR CAL, 

SATIEN  (4,3)  ,RVCAL,TVEL,GRCAL,NPOINT,SET,P  (3,3)  ,B  (3,3)  ,K1B,K1P, 
$RLlKE,KlD?H,K2EPH,K1DRKG, K2 DFKG, NPT  1 , XIC  (3, 1)  , BZERO, KlDELP, 
JK2JELP, K1DEIK, K2CELK, K  1  CP , K 2DF , IFLAG 
DIMENSION  HT  (3  ,3 )  ,FHT(3,3)  ,HPHT(3,3) 

01  CkuL  TRANS(H, HT,K1H,K2H,K1HT,K2HT) 

CAu L  MOLT (P,HT,FHT,K1F,K1F, KlHT,K2HT,K  1PHT,K2?HT) 

CAa.L  MJIT  (H,PHT,HPHT,K1H,K2H,  K  1 PHT,  K2P  fcT  ,  K1  H  PH  T,  K2HPHT) 

CALL  ADD(HFHT,E,B,K1HPHT,  K2H  PHT, R  IE,  K2E) 

RETORN 

ENj 

C 

C 

S  Ud^OUTINE  RINCVA  (XH  AT  ,  Z  MEA  S,ZH  AT,  R  N  0) 

C  SUBROUTINE  RINCVA  CALCULATES  THE  INNOVATION  FOE  THE  NLPI 

C  PROCESSOR. 

C  03 MON  F(3,  3)  ,G  (3,  1)  ,H  (3,3)  ,  R  (3 , 3)  ,Q  ( 1  ,  1 )  ,  K  IF  ,  K  1G  ,K2G  ,  K1 H  ,  K2n  ,  K  1 Q  , 
»K1i,RKGM  (3, 3)  , K 1FKIX , K2FHIX , K 1  PHI Z ,K 2PH1Z , 

JKlRKGM,  K2RKGM,RMASS,F11,F22,FKY,FKS1  ,F L E , ALPH A ,? Z ER 0 ,DT, NK ECR , 
iNR&CD,IPOSTA,IPOENC,IP 1ST A , IR 1END,I D AT S A ,IDATEN , N CAT  A, N R CAL, 

IATIEN  (4,3) ,RVCAL,TVEL,GFCAL,NPOIMT,SET,P  (3,3)  ,  B  (3 , 3)  ,  K  1c  ,K  IP, 
JRLIKE ,K1DPH,K2CRH,K1DRKG,K2DRKG,NPT1,XIC  (3, 1)  , BZERO, FIDEL?, 
JK2*)ELP,K1DEIK,K2DELK,KlCF,K2DF,  IFLAG 
DIMENSION  XHAT  (3,NPOINT) 

DIMENSION  ZMEAS  (3, NPOINT) 

DIMENSION  RNO  (3,NPOINT) 

DIMENSION  ZHAT  (3  ,  NPOIN  T) 

01  CA uL  MULT(H, XHAT, ZHAT, K 1 H, K2 H , 3 , NFOI NT ,K 1 ZHAT , K2 Z HAT) 

CALL  SOB  (ZMEAS  , ZHAT, PNU, 3, N POINT, K 1 R NU , K 2RNU) 

RETURN 

END 

C 

c 

SUBROUTINE  RM LP I  (R NU ,SC ,S D) 

C  SUBROUTINE  FKLFI  CALCULATES  THE  VALUE  OF  THE  LIKELIHOOD  FONC- 

C  TION. 

COMMON  F  (3 , 3 )  ,G  (3,1)  ,H  (3,3)  ,R(3,  3)  ,C  (1,  1)  ,  K1F,  K 1G  ,  K2G,  K1 H,  K2H  ,  K  IQ  , 
$K IF, RKGR ( 3, 3) , K 1PHIX , K2PH IX , K 1PHI Z, K2PHIZ , 

$K1MKGM,K2RKGM,PPASS, FI  1 , F22 , FRY , RK SI ,RLE , ALPH A, RZ ERO ,DT, NR ECR, 
$NRc,CD,  IPOSTA,IROEND,IP1STA,  IR  1  END,  ICAT  SA  ,  IDA  TEN,  N  DATA,  NRCAL, 
lATiiN  (4,3)  ,  RVC  AL,TV  EL,  G  RC  AL,  NFCI  NT,  SDT  ,  P  (3,  3)  ,B(3,3)  , RIB, KIP, 
$ELIKE,K 1DPH, K2CPH,K1DRKG,K2DRRG,NPT1,XIC  (3,1) , BZERO, KlDELP, 
jK2D£LP,K1DELK, K2DELK, K 1DF,K 2DF,I FLAG 
DIM  EN  SIGN  RNU  (3,NPOINT) 


n  n 
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DIMENSION  C  (3,  1) 

DIMENSION  ATC  (1,1) 

DIMENSION  AT  ( 1 ,3) 

DIMENSION  A  (3,  1) 

DIMENSION  BI  (3,3) 

DOUBLE  PRECISION  El 
DOUBLE  PRECISICN  NKAREA  (9) 

DOUBLE  PRECISION  BIDP(3,3) 

DOUBLE  PRECISICN  DETBD P 
DOUBLE  PRECISICN  E£P(3,3) 

DOUBLE  PRECISICN  VECTOR  (3) 

DOJBLE  PRECISICN  DID? 

DOUBLE  PRECISICN  D2DP 
DOUBLE  PRECISICN  BDP1(3,3) 

DOUBLE  PRECISION  FACTOR 
01  SUrt=0.0 
SC=0.0 
SD-0.0 
1DGT- 0 

CALL  SPTODP  (E  ,  K 1 E  ,  K 1 3  ,  ED P) 

PACTOR=1.00D  CO 

CALL  SCALDP  (9 DP , FACTOR , BDP1 , K IB , K 1 E , K 1 EDP  1,K 2BDP1) 

CALL  LINV1P(BDF,K1B,Kl3,BIDP, IDGT, HK AREA, IER ) 

CALL  DPTOSP  (9IDP  ,K1B,K1B,BI) 

IJJB=4 

D1DP=0.00D  00 

CALL  LINV3F(BDP1,VECTOR,IJOB,K1E,K  1B.D  ID  P,  D2D  P  ,  W  K  ARE  A  ,  IER ) 

IF (D2DP  .GE.  C.CCI  00)  GO  TO  50 

D2 JP--D2DP 

EX=2. 00D  00**  D2EF 

DErBDP=D1D?/EX 

GO  TO  60 

50  EX=2. 00D  00**D2DP 

DEr3DP=Dl DP*EX 
60  DE1'9=DET3DP 

S  1  -  ALOG  (DETB) 

H RITE  (6,6  5)  SI 

65  FORMAT (//,'  • ,12X, 'NATURAL  LOG  OF  THE  DETERMINANT  OF  B=»  , 

$  IX,  E15.8) 

DO  100  1=1 , NPC3NT 
DO  80  J=1 ,3 
A(J,1)=RNU (J,I) 

80  CONTINUE 

CAuL  MULT  (BI,A,C,3,3,3,  1,K1C,K2C) 

CA  i.L  TR  AN  S  (  A,  AT, 3,  1,KlAl,  K2AT) 

CALL  MOLT ( AT, C , ATC ,K 1 AT ,K 2AT, K 1C, K2C,K  1ATC,K  2 ATC) 

CONST=A  TC (1,1) 

SC=SC*CONST 
SD=SD*A  LOG (DETB) 

SUM=SUM+  (CONST  +  ALOG (DETB)  ) 

100  CONTINUE 

SD-  (“  1 .  0/2.  C)  *SE 
SCM-1. 0/2.0)  *SC 
RLIKE  = (“ 1 . 0/2. 0)*SUS 
RETURN 
END 


SUBROUTINE  PARPHI (EPMTX,DPUIX) 

C  SUBROUTINE  PARPHI  CALCULATES  THE  PAFTIAL  DERIVATIVE  OF  THE 
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C  STATE  TRANSITION  MATRIX  (PHI)  WITH  RESPECT  TO  THETA. 

COMMON  F(3,3)  ,G  (3,1)  ,H  (3,3)  ,R(3,3)  ,C(1#1)  ,K  1F,K  1S,K2G,K1H,  K2H,  Kly, 
$Kli,RKGH(3,  3)  , R 1PHIX, K2PHIX , K 1PHIZ , K2PHIZ , 

$K 1 H KGM , K2  RKGM, R HASS, F1 1 , P22, FKY , RK  SI , RLE ,  ALPH A, RZERO , DT, NR ECR, 

$N  R  ECD , IROSTA, IR  OE  N  E, IR 1ST A, IR 1  END, I  DAI SA , IDA  TEN, N  DATA, NRCAL, 

S ATTEN  (4,3) , RVCAL,TVEL,GRCAL, NPCI NT , SDT , P  (3, 3) , B (3,3) , K1B, K IP, 
$RLIKE,K1DPH,K2DPH,K1DRKG, K2 DR  KG , NPT 1 ,X IC  (3 , 1 ) , BZE FO, K1DELP, 

$K2 jELP , K1 DE1K, K2DELK , K 1  DP, K 2DF, I  FLAG 
DATA  ERROR/. OCCCOI/ 

DIMENSION  PRD  (3,3) 

DIMENSION  PRD  1  (3,  3) 

DIMENSION  TERM  (3,3) 

DIMENSION  SUM  (3,3) 

DIMENSION  DFHTX  (K1DF,K2EF) 

DIMENSION  DPHIX  (K  1  DPH  ,  K  2DPH) 

DIMENSION  PRD2  (3,3) 

DIMENSION  DIFF  (3,3) 

01  CA^L  SCALAR (F,C.O, DPHIX, KIP, KIP, K1DPH,K2DPH) 

CAi-L  SCALAR  (F ,  1 . 0  ,  PRD  ,  K  IP,  K  1 F,  K  1  PR D,  K2  £RD) 

DO  100  1=1, 20C 

CO~F=  (1/P ACT (I) ) *  (DT**I) *1 

C A mL  MULT (PPD, DFHTX, PRD2,K1PFE, K 2PR C , K 1 DF, K2 DP, K 1 PBD2,K2PRD2) 

C  A*.L  SCALAR  (PFD2,COEP,  TERM,  K1  PPD2,  K  2  EF  L2,  K  1TERK ,  K  2TESM  ) 

CAi.L  SCALAR  (DPHIX,  1 . 0  ,  S  UM  ,  K 1  DPH  ,  K2  DPH  ,  K 1  SUH,  K2  SUM ) 

C  Ai.L  ADD  (TERM, DPHIX, DFH IX ,K 1TEFH, K2TEPK, K1 DPH, K2DPH) 

CkuL  SUB(DPHIX,£UM,DIPP,K1DPH,K2DPH,K1EIPP,K2DIPP) 

C  ALL  MULT  (F , E FD , PFD 1 , K 1 F , K 1 P , K 1 PR D , K 2P F D, K 1PRD 1 ,  K 2P RD 1 ) 

C  A  i,L  SCALAR  (PRtl,  1.0,PRD,K1PFD1  ,K2PRD1  ,K  1PFD,  K2PR  D) 

DO  SO  J=1 ,K 1DIPF 
DO  50  K=1,K2BIFF 

IF  (AES  (DIFF  (J,K) )  -  ERR  OR )  50 , 50 , 1 00 
50  COST I NU  F 

GO  TO  200 
1G0  COS II NU  i 

WRITE  (6, 1  10) 

110  FORK  AT  (//, *  * ,  10X, *200  ITERATIONS.  NC  CO NV ERG £N CE  .  • ) 

200  RErUP  N 
END 
C 
C 

SUBROUTINE  PARF  (THETA, DFHTX) 

C  SUBROUTINE  EAPF  CALCULATES  THE  PAFT1AL  DERIVATIVE  OF  THE  SYSTEM 

C  MATRIX  WITH  RESPECT  TO  THETA. 

COMMON  F(3,3)  ,G  (3, 1)  ,H  (3,3)  ,R(3,3)  ,Q  (1,  1)  ,K1F,  K1G,K2G,  K1H,  R2H,  K1Q, 
SK 1 R, RKGK (3,3) , R 1PHIX , K 2PHIX , K 1PHIZ , K2P F IZ , 

SKlRKGM,K2RKGM,FHASS,F1 1 , F22 , FRY , RK SI , RLE , ALPH A , RZERO , DT, NR ECR , 
iNRiiCD,  IROSTA,  IROIND,  IR1STA,  IR  1  END,  I DATSA  ,1  DA  TEN ,  NDATA ,  NRCA  L, 

S  ATTEN  (4,3)  ,  PVC  AI,TV  EL,  G  FCAL  ,  N  POINT,  SDT,  P  (3,  3)  ,  B  (3, 3)  ,  K  IB ,  K  IP, 
JRL1KE,K1DPH,K2DPH,K1DRKG,K2DRKG,NPT1,XIC  (3,1) ,BZEFO, K1DELP, 

$K2 j£LP, K1 DELK,R2DELK,K1DF,K2DF,IFLAG 
DIMENSION  DFHTX  (K1DF,K2EF) 

DIMENSION  A  (3,3) 

A(  1,  1)  =  ((-2.0)  *P22)/(RHASS*T¥EL) 

A  ( 1 ,2)  =  (2.0*F22)/PHASS 
A  ( 1 , 3)  =0.  0 
A (2, 1) =0. 0 

A  (2,2)  =  (RKSI*TVEL)  /(2.C*(THETA**2.C)  *F  1 1  *  (PLE**2. 0 ) ) 

A  (2,3)  =0.  0 
A  (J,1)=0.0 
A (3, 2) =0. 0 


-  .)■>  « 'A 
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A (U ,3) =0. 0 

CALL  SCALAR  (A  ,  1.  0,  DFHT X , 3 , 3 , K IDF, K2 CF) 

RETURN 

END 


SUBROUTINE  PDP  (DFHTX , D ELTAP , l ELIA K) 

COMMON  F  ( 3f  3)  ,  G  (3, 1)  ,H(3,3)  ,R  (3,3)  ,C  (1  ,  1)  1 F, K  1G,  K2G,  K1  H,  K2B,  K1  a, 

$K1R, FROM (3 , 3) , K1PHIX,K2PHIX,K1PH1Z,K2PHIZ, 

SKlBKGM, K2HKGM ,RMASS,F  1 1 , F22 , R KI , R KSI ,R  LE,  AL  PH  A,  RZ  ERO  ,  DT,  NPLCR, 
SNRiSCD,IR0STA,IF0END,IR1STA,IR1END,  ID  ATS  A,  ID  ATEN,  ND  AT  A,  NR  CAL, 

SATTSN  (4,3)  , R VC A1 ,T VEL , GRC AL , NPOINT ,SET ,P  (3,3)  ,B(3,3)  ,K1B,K1P, 
*RLIKE,K1DPH,K2DPH,K1DRKG,K2DFKG,NPT1,XIC (3,1), BZEEO, K1D2LP, 
JK2jZLP,K1DELK,K2DELK,K1EP,K2DF,IFIAG 
DIMENSION  RI  (3,3) 

DIMENSION  HTR I  (3,3) 

DIMENSION  DELTAK  (3  ,3) 

DIMENSION  DFHTX  (3,3) 

DIMENSION  HT  (3  ,3) 

DIMENSION  DPHTXT  (3,3) 

D  1.1  CN SION  HP  (3,3) 

DIMENSION  RIH P  (3,3) 

DIMENSION  PRD  1  (3,3) 

DIIZNSION  FT  (3,3) 

DIMENSION  All  (3,3) 

DIMENSION  PRD2  (3,3) 

DIMENSION  PRD 3  (3, 3) 

DIMENSION  A21  (3,3) 

DIMENSION  RIH  (3,3) 

DIMENSION  HF.IH  (3,3) 

DIMENSION  PRD«  (3,  3) 

DIMENSION  A22  (3,3) 

DIMENSION  A  12  (3,3) 

DIMENSION  DELTA?  (3,3) 

DOUBLE  PRECISION  FHI(6,6) 

DIMENSION  ZDELP  (6,6) 

DOUBLE  PRECISION  RIDP(3,3) 

D0J5LE  PRECISION  RDP(3,3) 

DOUBLE  PPECISICN  W  K  A  R  E  A  ( 9 ) 

NAM  ELI ST/DI MEN/M 

C A i.L  TRANS(F,  FT,K1F,K1F,K1FT,K2FT) 

CAuL  SPTODP (R, FIR, KIR, REP) 

I D UI*0 

CALL  LINV  IF  (? E F ,  K  1 P ,  K 1  R ,  SI  D  P,  I DGT  , ¥K  AR  i  A  , I FR  ) 

CALL  DPTOSP  (RIDF,K1E,K  IF,  FI) 

CALL  TRANS(H,HT,K1H,K2H,K1HT,K2HT) 

C  Ai<L  1RANS(DFHTX,  DFHT  XT,K1F,K1F,K  1EFH,K2DFH) 

CA  a,L  HOLT  (H,P,HP,K1H,K2H,K1  F,  K 1  P  ,K  1  HP,  K2HP ) 

CALL  MULT  (R I , HP ,R I  HP , K 1R, K 1R , K 1HP, K2HP ,K 1 RIH P, K2 BIHP) 

CA^L  HOLT (HT,RIHF,PFD1,K1HT,K2HT, K IF 1HP, K2R1HP, KlPRD  1  ,K2PRD 1) 

CAuL  SUB (PR D1 , FT, A 1 1 , K 1PRD 1 , K2PRD 1 , K  1 A  1 1 , K 2A 1 1) 

CA uL  HULT (P,DFFTXT,PRD2,K1P,K 1P,K1DFH, K2DFH, K 1 PP D 2, K 2PRD2) 

CALL  HULT  (DFHTX ,P ,PRD3 , RI DFM, K2EFB , KIP ,K 1P,K 1 PRD3, K2PRD3) 

CAuL  A DD (PRC2 , FPD3 ,A21,K1PRD2,K2PRD2,K1A21,K2A21) 

CA*.L  HULT (R I, H,RIH, KIR , K1R ,K 1 H,K2H, Kl RIH, K2RIH) 

CALL  HULT  (HT,  RIH,HRIH,  K  1HT,  K2HT,  K  1RI H,  K2RIH,  K  tHPI  H,  K2HRIH) 

C  Ai.L  HULT  (P,  HRI H,  PRD4 ,  K  IP, K  IP,  K 1HRIH , K  2HRIH,  K1 PRD4,  K2PRD4) 

CALL  SUB(F,PRD4,A22,K1F,K1F,K1A22,K2A22) 

CA-L  SCALAR (F, 0 . 0, A12 , K 1F,K IF, K1A12 , K2 AI 2) 

DO  20  1=1  ,3 


no  no 


00  20  J=1  ,3 
ZDELP  (I ,  J)  =  A1 1  (I,J) 

20  CONTINUE 

DO  40  1=1,3 
DO  40  J=1 ,3 
M  =  I 

N= J+K2A 1 1 

ZDJLP  (M,N)=A12  (I,J) 

40  CONTINUE 

DO  60  1=1,3 
DO  60  J  =  1  ,3 
H=I  +  K  1 A  1 1 
N=J 

ZDiiLP  (M,  N)  =A2  1  (I,J) 

60  CONTINUE 

DO  80  1=1,3 
DO  80  J*1,3 
M=I+K 1 All 
N=J*K2A  Vi 

ZDcLP  (9 , N ) =A22  (I, J) 

80  CO&TINUE 

WRITE  (6,100) 

100  FORMAT (//, •  ' , 12X, ' Z  MATRIX  FOP  SOLUTION  TO  PARTIAL (P) /• ff 
$• PARTIAL (THETA) •) 

CAkL  OUTPUT (ZDEIP,H,M) 

IRaTE  (6,DIH2N) 

CALL  STRZ (ZBEIP,ET,PHI,M, KlPHI,K2PHI) 

WRITE  (6,120) 

120  FORMAT (//, '  ',12X,'STK  FOR  SOLUTION  TO  P APTIA L (P) /• , 

S' PARTIAL (THETA) ') 

CALL  OUTDP(PHI,K1PhI,K2PHl) 

CALL  RICATI (PHI ,K 1  PHI , DELTA  P,K1DELP) 

CALL  TRANS (H, HT ,K 1H, K2 H,K 1HT, K2HT) 

CALL  HOLT  (HT ,  R I , HTR1 , K 1 HT , K2H T, F  1 R, K  1R, K  1HTP.I  ,  K  2HTRI ) 

CALL  MULT (DELT AP, HTP I , DELTA K , K 1 DELP , K2 1 ELP ,K 1 HT RI , K2 HT RI , K1 DELK , 
$K2  JBLK) 

RE1UFK 

END 


SUBROUTINE  XSEN (PHIK, K 1  PH IK , DFBTX, LI LT  AK , XHAT, ZF EAS, DXHAT) 
SUBPOUTINE  XSEN  CALCUALTES  THE  PARTIAL  DERIVATIVE  OF  XH AT 
WITH  RESPECT  TO  THETA (TH  E  VECTCF  CF  UNKNOWNS. 

COMMON  F  ( 3,  3)  ,6(3,1)  ,H{3,3)  ,  R  (3 , 3)  ,Q  (1,1)  ,K1F,K1G,K2G,K1H,K2H,K1Q, 
$K1R,RKGR (3,3) , K 1 PHIX , K 2  EH  IX , K 1PHI Z , K 2PHI Z, 

SKlRKGK, K2RKGM,RMASS,F1 1 ,F 22 , P K Y , RKSI ,RI E, A LPHA , PZ FRO, ET.  NRECR, 
$NGaCD,IR0STA,IF0END,IF1STA,IR1END,  ID  ATS  A,  ID  ATE  N,N  DAT  A,  NR  CAL, 

SATEEN  (4, 3)  ,  RVCAL,TVEL,  GRC  AL,  NPOINT  ,  SIT  ,  P  (3,3)  ,B(3,3)  ,K1B,K1P, 

SR LIKE ,K1DPH,K2DPHrK1DFKG,K2DRKG,NPT1,XIC (3,1), BZEFO, K1DELP, 
SK2DELP, K1DELK, A 2CEIK, K  1 EF , K 2DF , IFI AG 
DOUBLE  PFECISIC  N  FHIK(3,3) 

DOUBLE  PR  EC  I  SI  C  N  PHIKI(3,3) 

DOUBLE  PRECISION  BKAREA  (9) 

DOUBLE  PRECISION  DPKH  (3,3) 

DOUBLE  PRECISION  CFKHZ(3,3) 

DIMENSION  SPHIKI  (3,3) 

DIMENSION  DFHTX  (3,3) 

DIMENSION  DELTAS  (3,3) 

DIMENSION  XHAI  (3, NPOINT) 

DIMENSION  ZHEAS  (3, NPOINT) 
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DIMENSION  DXHAT  (3,NP0INT) 

DIMENSION  RB  (3 , 3) 

DIdEN SION  PBD 1  (3,3) 

DIMENSION  X1NPT  (3,3) 

DIdEN SION  SPHIK  (3,3) 

DIMENSION  RRGH  (3,3) 

DIdENSI ON  FKH  (3,3) 

DIdEN  SION  PBD  2  (3,3) 

DIMENSION  PRD3  (3, 3) 

DIdEN SION  PRD4  (3,3) 

DIMENSION  PPD5  (3,3) 

DIdEN SION  PBD6  (3,3) 

DIMENSION  PKHI  (3,3) 

DIdENSI ON  VXH  AT  (3,1) 

Did  ENSION  VDXH AT  (3,1) 

DIMENSION  VZ  (3,1) 

DIMENSION  TERM  1  (3,  1) 

DIMENSION  TEFN2  (3,1) 

DIMENSION  TERM3  (3,1) 

DIMENSION  SUM1  (3,1) 

DIMENSION  VXSEN  (3,1) 

DIMENSION  GAMMA  (3 ,3) 

DIMENSION  LAMBDA  (3,3) 

CA*.L  MOLT  (DELI  AK  ,H  ,  PPD  1,K  1DELK,  K2DELK,  K  1  fl,  K2U,  K  1PRD 1 ,  K2PED 1) 

C A i.L  SUB(DFMTX,FRD1,XIXFI,K1?,K1F,K1XIKP,K2XINP) 

CAi*L  DPTOSP (PAIR, 3, 3, SPHIK) 

PA  C 10  R=  0 . 0 

C  A»L  SCALAR  (F,PACTOP ,FM,K1F,K1F,K1RP,K2PM) 

DO  40  1=1 ,K1 F 
FB  (I,I)=1.0 
4  0  CONTINUE 

IDGT=Q 

CA-L  LINV 1 F (PHI K, K 1 PHIK, K 1PHIK, PHIKI,1DG1, HK ARE A, IE?) 

CAuL  DPTOSP  (PHIKI,K1PBIK,K1PHIK,SPHIKI) 

CALL  MOLT  (RKGR,H,RKGH,  K  1 R  KG  M,  K2FK  GB,  K  1  b,  K  2H,  K  1  RKGH,  K2RKGH) 

CA„L  SUB(  F,  PKGH,PKH,  KIP,  KIP,  K1FKH,K2FKH) 

CAi.L  SPTODP  (PKH,KlFKH,K2FKH,DFKH) 

IDGT=  0 

CALL  LINV 1F(DFKfi,KlFKH, K2FKH, DFKHI, I DGT, HK  AREA, IE R) 

C A„L  S0E(RR,SPHIKI,PRD2,K1RK,K2RB,K1FRD2,K2PRD2) 

CA„L  MOLT (PPD2,XINPT,FHD3,K1PFD2,K2PFD2,K1XINP,K  2XIN  P,K1PRD3, 
SK2PBD3) 

CALL  DPTOSP (DFKII,K1FKB,K2FKH,PKHI) 

CA„L  MULT(FKHI,PPD3,PRD4,K1FKH,  K2F KH  ,  K  1 PRD3,  K2 PRD3,  K1PRD4,  K2PRD4) 
CA«L  B'JLT  (SPHIK,PRD4,GAKMA,K1FHIK,K1PHIK,K1PFD4,K  2PRD4,  K 1  GAMA, 
$K2JAM  A) 

CALL  MOLT (PRD2 , DELTA K, PFD5, K 1PRD2, K 2FPD2, K1DELK, K2DELK,K 1PRD5, 
SK2PBD5) 

CA-L  MULT (FKHI,PRD5,PFD6,K1FKH,R2FKH,K1PFD5, K2PRD5,K  1PRD6,K2PRD6) 
CALL  H0LT(SPHlK,PRD6,lAP3DA,K1PHIK,K1PtIK,R1  PRD6,  K2PRD6, 
$KUABE,K2LAMB) 

DO  200  J=  1,  3 
VDXHAT (J,1) =0.0 
200  CONTINUE 

DO  b 00  1=1, NPC 1NT 
DO  220  J=  1,  3 
VXiAT  (J,  1  )  =  XU AT  (J,I) 

VZ  (J,  1)  *ZBEAS  (J,I) 

220  CONTINUE 

CAi,t  MULT  (SFBIK,VDXHAT,TEFM1,K1PRIK,K1FHIK,3,  1,  K  ITER  1 , R2TEB  1) 
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CALL  HOLT  (G  AHH  A  #¥XHAT  f  TEFM2,  K  1G  AH  A,  K2G  AH  A  ,3,  1  ,K  1TER2,  K2TER2) 

CALL  HOLT (LAHBDA,VZ,TEBH3 ,K 1 LAHB, K2LAHE , 3, 1 , K 1TEH 3,K2TER3) 

CALL  ADD  (TERN  1 ,T£RH2,SU Hi , K ITER  1 ,K2TER  1,K 1SUN  1.K2S0H  1) 

CALL  ADD(SUHl,TEPf!3,¥XSfN,KlSUHl,K2S0Hl,  R 1 VXSN  ,  K2VXS  N) 

DO  400  J=1,3 

DXtiAT  IJ,I)=VXSEN  (J,1) 

400  CONTINUE 

CAkL  SCALAR (VXSEN,1.0,VEXHAT,K1VXSN, K2VXSN, K1 VDHT.K2VDHT) 

500  CONTINUE 
RETURN 
END 
C 

c 

SUBROUTINE  GRACE  (EXHAT,  RNU, DELTA? , CELT AJ) 

C  SUBROUTINE  GRADE  COMPUTES  THE  GRADIENT  OF  THE  LIKELIHOOD 

C  FUNCTION. 

COMMON  F(3,3)  ,G  (3,  1)  ,H  <3, 3)  ,R  (3,3)  ,Q  (1,1)  , K  1 F,  R  1G,  K2 G,  K1  fl,  K2H,  K1  Q , 
$Kla,RKGM  (3,3) , K 1 PHIX , K 2 EHIX , K 1  PHIZ , K2PKIZ, 

SK 1nKGM,K2RKGM ,RMASS,F1 1 ,F22 ,RKY, RKSI  ,R IE, ALPHA, RZERO,DT, NRECR, 

SNR LCD, I  POST A, I R0END,IR1STA,IR1END, ID AT S A, ID ATEN, ND AT  A, NR CAL, 

SATTEN  (4,3) ,R VC AI,T VEL , GRC AL , NPOINT  ,SDT ,P  (3 ,3 )  ,B  (3 , 3)  , K1 B, K 1 P, 
SRLIKE,K1DPH,K2DPH,K1DRKG,K2DRKG,NPT  1, XIC  (3, 1)  , BZERO, K1DELP, 
SK2JELP,K1DELK, K2DPLK,K1DF,K2DP,IFLAG 
DIMENSION  DXHAT  (3, NPOINT) 

DIMENSION  DELTAP  (K1EEL?,K2DELP) 

DIMENSION  RNU  (3, NPOINT) 

DIMENSION  A  1  ( 3 , 1) 

DIMENSION  AIT  (1,3) 

DIMENSION  BI (3,3) 

DIMENSION  HI  (3,3) 

DIMENSION  A2  (3  ,3) 

DIMENSION  A3  (3,3) 

DIMENSION  DELTAE  (3 ,3) 

DIMENSION  A 4  (3,3) 

DIMENSION  A5  (3 ,3) 

DIMENSION  A6  (3,  1) 

DIMENSION  HA6  (3,1) 

DIMENSION  A7  (3 , 1) 

DIMENSION  A8  (3 , 1) 

DIMENSION  TERM  1  (1,1) 

DIMENSION  TERN 2  (1,1) 

DOUBLE  PRECISION  WKARE A  (9) 

DOUBLE  PRECISION  RDP(3,3) 

DOUBLE  PRECISION  BIDP{3,3) 

DIMENSION  DELT  A V  (3,1) 

01  SUS=0.0 

IDGT=  0 

CALL  SPTODP  (B,K1B,K1B,BDP) 

CALL  LINV1F(BDF,K1B,K1B , BID P, IDGT ,BKARIA,IER) 

CALL  DPTOSP  (BILP,K1B,K  IE, BI) 

CALL  THANS(H,HI,K1H,K2H,K1dT,K2HT) 

CALL  HULT  (DELTA  P,  HT,  A2 ,  X  IDE LP,  K2DELP, K  1HT,K2HT,K  1  A2,K2A2) 

CALL  HULT (H, A  2, CELTAB , K 1 H,K 2 H ,K 1 A2 ,K2A2, K1 DEL3,K2D8LB) 

CALL  HULT  (DELT AE ,BI, A3 , K 1DBLB, K 2DELE ,K IB, K 1 B, K 1A3,K2 A3) 

CALL  HULT  (BI ,A3,A4,K1B,K1B,K1A3,K2A3,K1A4,K2A4) 

CALL  HULT  (BI,DE1TAB, A 5 , K1B, R IB, K 1  DELE, K2DELB, K 1 A5 ,K2 A5) 

CAi.L  TRAC  E(A5,T3,K1A5) 

DO  100  I«1, NPOINT 
DO  20  J«1,3 
A 1  (J,  1)  «HNO  (J,I) 
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A6(J,  1)  *DXHAT  (J,I) 

20  CONTINUE 

CALL  TRANS  (A1,A1T,3,1,K1A1T,K2A1T) 

CALL  HOLT  (H,A6,HA6,K1H,K2H,3,  1 ,  K  1H  A6  ,K2H  A6) 

CALL  NEG (HA6,DELTAV,K1HA6,K2HA6,K1 CELT ,  K2DEL  V) 

CAuL  HOLT  (BI,DELTAV,A7,F1B,K1B,  K  UE1V,  K2DELV,K  1  A7,K2A7) 

CALL  HOLT  (A1 T,  A7,  TEBH 1 ,  K1 A1  T,  K2A IT,  K  U7,  K2A7  ,  K  1TER1  , R2TBR 1 ) 

CALL  HOLT  (A4,  A  1 ,  A8  ,K  1 A  4  ,K2A4 , 3 , 1  ,K  1 A8,  F2  A8) 

CALL  HOLT  (A  IT  ,  A8,TER  32  ,  K1  A  IT,  K2  AM,  K  Ue,  K2A8, K1TER2,  K2TEB2) 

T  1=TE  BH 1  (1,1) 

T2-TERM2  (1,1) 

SUfl=SOM*Tl*  ( (1. C/2.0)*  (13-T2)  ) 

100  CONTI NUE 

DELTAJ-SOH 

flEXOBN 

END 

C 

C 

SUdROUTINB  FISBS  B  (FNU ,DXHAT, DELTAP,  EELJ2) 

C  SUBROUTINE  FISHER  COMPOTES  THE  FISHER  INFORMATION  MATRIX. 

C  THE  FISHER  IKFC5MATICN  MATRIX  IS  THE  SECOND  PARTIAL  DERIVATIVE 

C  OF  THE  LIKELIHOOD  FUNCTION  HITH  RESPECT  TO  THETA. 

COMMON  F(3,3)  ,G  (3,  1)  ,H  (3,3)  ,R  (3,3)  ,C  (1,  1)  ,K  IF, K 1G, K2G, K1 H, K2H, K 12 , 
»K IE, RKGH (3, 3) , K 1PHI X, F 2  PHI i ,K 1  PHI Z , K2P HI Z, 

SK  IdKGM,  K2RKGM  ,RH.»?S,F  1 1,F22,  PKY,  RK5I  , RLE,  ALPHA,  RZERO,  DT,  NRECR, 
$NRECD,IPOSTA,irOEND,IPlSTA, IR 1 EN D, It AT SA ,IDAT2N , NDAT A , NR CAL, 

RATTEN  (4,3)  ,  R  VC  A  WTVEL,  GRC  A  L  ,  NPOI  NT  ,SET,P  (3,3)  ,  B(3,3)  ,K1B,K1P, 
SfiLIKE,K1DPH,K20i'H,KlwPKG,K2DRKG,NPTl  ,XIC  (3,1)  ,  BZERO,  K1DELP  , 
SK2JELP, K 1DELK , K2 DELK , K 1 DF , K  2DF, IF! AG 
DIMENSION  RNU  (.*„  <*  OINT) 

DIMENSION  D  '*  -.1  AT  ( 3  ,  NP  0 1  NT ) 

DIMENSION  DELTAP  (K 1 D ELP , K 2D  EL P) 

DIMENSION  DELTiE(3,3) 

DIMENSION  :5  N  £  3  (3,3) 

DIMENSION  HI  (3,3) 

DIMENSION  DELN 0  (3,1) 

DIMENSION  DELNUT  (1  ,3) 

DIMENSION  DXHAT1  (3,1) 

DIMENSION  RNO  1  (3,1) 

DIMENSION  RN1T  (1,3) 

DIMENSION  TERM  1  (1,1) 

DIMENSION  TEPM2  (1,  1) 

DIMENSION  TERM 3  (1  ,  1) 

DIMENSION  A1  (3,3) 

DIMENSION  PRD  1 1  (3,1) 

DIMENSION  PRD2 1  (3,1) 

DIMENSION  PRD22  (3,1) 

DIMENSION  PPD23  (3 ,  1) 

DIMENSION  PRD31  (3,1) 

DIMENSION  PRD32  (3,1) 

DIMENSION  PRD33  (3,1) 

DIMENSION  PRD41  (3,3) 

DIMENSION  PPD42  (3,3) 


DIMENSION  PRD43  (3,3) 
DIMENSION  BI  (3,3) 

DOUBLE  PRECISION  VKAREA (9) 
DOUBLE  PRECISION  BDP(3,3) 
DOUBLE  PRECISION  SIDP(3,3) 
SUM* 0.0 
IDGT*0 


CALL  SPTODP  (E, K1B,K1B,BDP) 

CALL  LI  MV  IF (BEF ,K 1 B, Kl B ,BIDP, IDGT, R KAP IA , IBB ) 

C Ai.L  DPTOSP  (BILF,  K1B,F1B,BI) 

CALL  IRANS(H,H1, K  1H,  K2H,K  1HT,  K2HT) 

CALL  NEG  (H,HNE€,K1H,K2H,K1HNEG, K2HNEG) 

CAi.L  BOLT  (DELTAF,HT,A1,K1  DELP  ,K2DELP,K1HT,K2HT,K1A1,  K2A 1 ) 

CALL  BOLT  (H, A 1 #DELTAB , K1H, K2H ,K1 A 1, K2 A  1, K 1DEL B# K2DELB) 

DO  100  1=1,  NPCINT 

DO  20  J=1 ,3 

DXHAT  1  ( J,  1)  =DXH AT  (J,l) 

EN  J1  (J,  1)  =FN 0  (J,I) 

20  CONTINUE 

C  COMPUTE  THE  1ST  TERM  OF  THE  INFORMATION  MATRIX 

CAi.L  BULT  (HNEG,EXHAT1,DELNU,K1HNEG,K2HNEG,3,  1,  KID  NO,  K2DNU) 

CA-L  BULT  (BI,DELNU,PBD1 1 , K 1 B , K 1 E , 3 , 1  ,K1PD1 1,  K2PD11) 

CALL  TRANS ( DEL NU , DELNO T , K 1 D NO , K 2D NO, K 1BN OT, K2DNUT) 

CALL  BOLT (DELN OT, PRD1 1 , TERM  1 ,  K 1DNUT , K2 E NUT, K 1  PD1 1 , K2 PD1 1 ,  K 1  TER  1 
SK2TE.R  1 ) 

C  COMPUTE  THE  2NE  TERM  OF  THE  INFORMATION  MATRIX 

CAuL  MOLT  (BI,DELNU ,PFD 2 1, K 1 E, K IB, K ID NU, K 2DN0, K 1PD21 , K2PD2 1 ) 

CALL  HULT(DEL1AE,PRD21 ,P3D22, K1DELB, K2EFLB,K1 PD21 ,K2  PD21,K1P L22 
SK2PD22) 

CAi.L  BULT(Bl,PRt22,PBD23,K1E,KlE,KlPE22,K2PD22,K1  PD23,  K2PD23) 
CALL  TRANS(ENU1,BN1T,3, 1, K 1 PN IT, K2RN  IT) 

C Aj.L  MULT(RN1T,FRD23,T£RM2,  K 1 FN IT , K2RN IT, Kl PD23 , K2PD23 , Kl TEF2 , 
SK2TER2) 

C  COMPUTE  THE  3R C  TERM  OF  THE  INFCFHA1ION  MATRIX 

CALL  MOLT (BI, DELNO , PFD 3 1 , K 1 B, K IB, K ID  NO, K 2DNU, K 1PD3 1 , K2PD3 1 ) 

C  A*.L  BOLT  (DEL1AE,PFC31,PRD32,K1EFIB,K2L£LB,K1  PD3 1  ,K2  PD3 1,K1PL32 
SK2PD32) 

CALL  MULT (B I , P B E32 ,PR D 33 , K 1 E , K 1 E, K 1 PE3 2 , K2PD3 2 , Kl PD33, K2PD33) 
CALL  MOLT (R Nl T, FRD33 , TE R13, K 1 PN IT, K 2FN  1T,K  1PD3 3, K  2PD33  ,K  ITER  3, 
$K2  IER  3) 

C  COMPUTE  THE  4TH  TERM  OF  THE  INFORMATION  MATRIX 

C  Ai.L  BOLT  (El,  DEITAB,PRD41  ,K  IE , K 1 B, Kl EE  IB , K2DELB,  K 1 PD4  1,K2PD41) 
CALL  BOLT (DELTAB,PRD4 1, FRD4 2, K 1DELB, K2DEL E, K 1PDU 1 ,K2PD41 , K  1PD42 
SK2PD42) 

CALL  MOLT (BI, PBD42 , PPD4 3, K 1 B, K IB, K 1 PD42, K2PD4 2, K 1PD4 3, K2PD43) 
CALL  TRACB(PRD43,T4,K1PE43) 

I1=TERM1  (1,1) 

T2=TEPM  2(1,  1) 

T3=TERM3  (1,1) 

PSJM=T1-T2-T3-  ((1.0/2.C)*T4) 

SUM  =  SUM  *P  SUM 
ICO  CONTINUE 
DE*,J2=S0M 
RETURN 
END 
C 
C 

SUBROUTINE  STFF (TH ETA, DELTA J , DELJ2 , ETHET A,TH ETAN) 

C  SUBROUTINE  STEP  CALCULATES  THE  NEW  VALUE  OF  THETA  (THETAN) . 

DTdETA*  (“1.0/BELJ2)*DEL1AJ 
TH3TAN=DTHETA*TBETA 
RETURN 
ENJ 
C 
C 

SUBROUTINE  BULT  (A , B,  D,  K  1A ,K 2 A , K IB , K2E,  M  D, K2  D) 

DIMENSION  A  (K  1 A  ,E2 A)  ,  B  (K  IB,  K2B)  , D  ( K  1 1,  K2B) 
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IF  ( K2  A  .NE.  K 1 B) GC  TO  40 
DO  20  1  =  1, K1A 
DO  20  J=1,K2B 
L=J 

DO  20  K=  1 ,K2A 
L=i.*1 

IF(L  .HE.  1)  GO  TO  20 
D(I,J)  =0 

20  D(I,J)=D{I,  J)  ♦  (A(I,K)*3  <X,J)) 

KlJ=K1A 
K2D=K2B 
GO  TO  80 
40  WRi.TE  (6,60) 

60  FOR!1AT(1HO,  'MA1RICES  A  AMD  B  ARE  NOI  CONFORMABLE* ) 
80  RETURN 
ENJ 


SUBROUTINE  SCALAR  (A, B,C,K1A,K2A,K1C,R2C) 
DIMENSION  A  (K  1  A  ,K2A)  ,C  ( K 1  A, K2A) 

DO  10  1=1 ,K1A 
DO  10  J=1,K2A 
10  C (I,J) =B*A (I, J) 

K 1 J=K  1 A 
K2C= K2 A 
RETURN 
INJ 


SUBROUTINE  TRANS  (A,B,K1A,K2A,1<  IE,  K2E) 
DI.1  EKSION  A  (K1A,K2A)  ,B  (K2A,K1A) 

DO  20  1=  1 , K  1 A 
DO  20  J=1 ,K2A 
20  3  (J,I)  =A(I,  J) 

K 1d=K2A 
K2B  =  K  1A 
RETURN 
ENJ 


SUoPOUTINE  NE3  ( A ,  B,  K  1 A ,  K2 A,  K  1 B,  K2P) 
DIMENSION  A  (K1A,K2A)  ,B  (K1  A,K2A) 

DO  10  1=1, K1A 
DO  10  J=1, K2A 
10  B (I,J) =0- A  (I,J) 

K1b=K1A 

K23=K2A 

RETURN 

ENJ 


SUBROUTINE  ADC  (A,B,C,K1A,K2A,K1C,K2C) 
DIMENSION  A  (K  1  A,K2A)  , B  (K 1 A ,  K2A)  ,C  (K  1  A,K 2  A) 
DO  20  1=1,  IMA 
DO  20  J=1,K2A 
20  C  (I,J)=A  (I,J)  4B  (I,J) 

K 1C=X  1 A 
K2C=K2A 
RETURN 
ENJ 
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c 
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SUBROUTINE  SUE  (A,  B,C  ,  K  1A,  K2  A,  K  1C,  K2C) 
DIMENSION  A  (K  1  Af  A2A)  ,B(K1A,K2A)  ,C(K1A,F2A) 
DO  20  1*1, K1A 
DO  20  J*1 , K2A 
20  C(I,J)=A(I,J)-E  (I,J) 

K1C*K1A 
K 2C=K 2A 
RETURN 
END 
C 
C 

FUNCTION  FACT (K) 

FACT* 1.0 
DO  20  1*1  ,  K 
RI*I 

20  FACT=FACT*EI 
PETURN 
END 

c 

c 

S UD SOUTINE  TRACE  (A,B,  K1A) 

DIMENSION  A(K1A,K1A) 

sua=o.  o 
DO  20  1*1, K1A 
SU.1  =  SUM*A  (1,1) 

2C  CONTINUE 
B=3UM 
RETURN 
END 
C 
C 

SUBROUTINE  0UTPUT(A,K1 A,K2A) 

DIMENSION  A  (K  1 A  ,K2A) 

DO  20  1*1, K1A 

WRITE  (6,30)  (A  (I,J)  ,J*1,K2A) 

20  CONTINUE 

30  FORMAT ( •  ' , 12X, t  |2X,£1 5. 8) ) 

RETURN 

END 

C 

C 

SUBROUTINE  OVER  (RM ATR X , CdK , K1 ,K2,IFLIGO) 
DIMENSION  EHATFX  (R1.K2) 

DO  20  1=1 ,K1 
DO  20  J=1,K2 

IF  ( RNAT  PX  (I ,  J)  .GE.  CHK)IFLAGO*1 
20  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  SFTCC  F  (X,K1X,K2X,Y) 

DIMENSION  X (K  1X,K2X) 

DOUBLE  PRECISION  I(K1X,R2X) 

DO  20  1*1, KlX 
DO  20  J*1,K2X 
Y  U,J)  *DBLE  (X  (I,J)) 

20  CONTINUE 
RETURN 
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END 


SUBROUTINE  DPTOSP  |X,K1X,K2X,Y) 
DOUBLt  PREC1SICK  X(K1X,K2X) 
DIMENSION  Y(K1X,K2X) 

DO  20  1*1 ,K1X 
DO  20  J=1,K2X 
Y  (X,  J)  *X  (I,  J) 

20  CONTINUE 
RETURN 
END 


SUBROUTINE  SCALDP (A,B,C,KlA,K2A,KlC,K2C) 
DOJBLE  PRECISION  A  (K  1A  ,  K2A)  ,  C  (K  1A  ,K2  A) 
DOJ3LE  PRECISION  B 
DO  10  1*1, K1A 
DO  10  J=1 ,  K2  A 
10  C(X,  J)  =  B*A(I,J) 

K  1 J  =K  1 A 
K2C=K  2A 
RETURN 
END 


SU3ROUTIN  E  TR AN SD  ( A , B , K 1 A , K2 A , K  1 B, K2 E) 
DOUBLE  PRECISION  A ( K 1  A ,K 2A>  , E  ( K2 A, K  1  A) 
DO  20  1*1, K1A 
DO  20  J=1 , K2  A 
20  B  ( J  ,  I )  =  A  ( I,  J) 

K 1 j*  K2A 
K23=K  1A 
RETURN 
ENJ 


SUBROUTINE  NEGCP  (A,B,K1A,K2A,  K1E,  K2 E) 
DOJBLE  PRECISION  A  (K  1  A  ,  K2  A)  ,  B  (K 1  A ,  K2  A) 
DO  10  1=1, K1A 
DO  10  J=1 , K2A 

10  B  (x,J)=DBLE  (0.  0)-A  (1,J) 

K1J  =  K1A 
K  2i=K  2A 
RETURN 
ENJ 


SUBROUTINE  ADEC P  (A,B,C,K1A,K2A,K1C,K2C) 

DOJ3LE  PRECISION  A  (K  1 A  ,  K2A)  ,B  (K  1 A  ,K  2  A)  ,C  (K  1  A,  K  2A) 
DO  20  1*1, K1A 
DO  20  J* 1 ,K 2A 
20  C  (I,J)*A  (I,J)  *E  (I,J) 

K 1C*K  1 A 
K2C*K2A 
REX0FN 
ENJ 
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SUBROUTINE  MULTDP  (A ,B  ,D  ,  X1A ,  K 2A,  K  1 B,  K2  E,  K  ID,  K 2D) 

DOJBLE  PRECISION  A  (K  1 A  ,  K2  A)  ,  B  (K1  B  ,  K2  E)  ,  D  (K 1 A  ,  K2B) 

IF  (K2A  .NB.  K 1  E)  GC  TO  40 

DO  20  1*1, KlA 

DO  20  J=1 ,K2B 

L=0 

DO  20  K= 1 ,K  2 A 
L=L*1 

IF(L  .NE.  1) GO  TO  20 
D  (I,  J)  =DBLE  (0.0) 

20  D(I,J)»D(I,J)*(»(I.K)*B(K,J)) 

KlD=K1 A 
K20=K2B 
GO  TO  80 
40  BRITE  (6,60) 

60  POBMAT(1HO, 'MATRICES  A  AND  B  ARE  NOT  CCNPORMABLE' ) 
80  RETURN 
END 
C 
C 

SUBROUTINE  SUBDP  (A ,S , C , K 1  A, K2 A , K 1 C  ,  K2C ) 

D0J3LE  PRECISION  A  (K  1 A  ,  K2A )  ,  B  (K  1 A  ,  K  2  A)  ,C  (K  1  A,  K  2  A) 
DO  20  1=1, KlA 
DO  20  J=  1 , K2 A 
20  C  (I,J)=A  (I,  J)  “E  (I  ,J) 

K  1C=K  1  A 
K2C=K2A 
RETURN 
END 
C 

c 

SUBROUTINE  TR ACED  (A, B,K  1A) 

DOJBLE  PRECISION  A(K1A,K1A) 

DOJELE  PRECISION  SUM,B 
SU.i=D3LE  ( 0.  0) 

DO  20  1=1  ,  K 1 A 
SU3=SUM+A  (1,1) 

20  CONTINUE 
B=3UB 
RETURN 
END 
C 

c 

c 

c 

SUBROUTINE  OUTD3L  (A, K 1  A ,K2 A) 

DOJBLE  PRECISION  A  (KlA, F 2 A ) 

DO  20  1=1, K1A 
DO  20  J=1,K2A 
IRiTE  (6,30)  A  (I,J) 

20  CONTINUE 
30  FOBHATC  *,12X,D15.8) 

RETURN 

END 

C 

c 

SUBROUTINE  0?  E  RDP  (RMAT  RX,CHK,  K 1 ,  K2, 1  FL  AGO) 

DOUBLE  PRECISION  RMATFX  (K 1 ,  K2) 

DOJBLE  PRECISION  CHK 
DO  20  1=1, K1 
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I 

I 

I 


DO  20  J=1 ,K2 

I F  (EH  ATRX  (I ,  J)  .GE.  CHK)IFLAG0=1 
20  CONTINUE 
RETURN 
£N  u 
C 

c 

DOUBLE  PRECISION  FUNCTION  PACTDP(K) 

DOJBLE  PRECISION  DREALI 
F ACTDP=DBLE  (1.0) 

DO  20  1=1, K 
RE  ALI=I 

DREALI=DBLE  (REALI) 

20  PACTDP=FACTDP*CRE All 
RETURN 
ENJ 
C 

c 

SUBROUTINE  OUTEP  (A  ,K  1  A  ,  K2A) 

DOJBLE  PRECISION  A  (K  1  A,  K2  A) 

DO  20  1=1, K1A 

WRITE  (6  ,3 0)  (A  (I,J)  ,J=1,  K2A) 

20  CONTINUE 

30  POaflATC  '  , 12  X  ,6  (2X,  D  15.8)) 

RETURN 

ENJ 

C 

c 

SUBROUTINE  PLOT  (X ,  Y ,  IPO  1  NT) 

DIMENSION  X (IPCINT) ,Y (IFOINT) 

WRITE  (6,  20) 

2C  FOaM  AT (/, •  ' , 4  OX , *  STATE  ESTIMATES') 

CA^L  WPLOT1 (X, Y,IPOINT,4,XHAT) 

RETURN 

ENJ 

C 

C 

SUBROUTINE  VS  N  E X  (X, K  X ,  R  MI N ,  FR  A X > 

DIMENSION  X  (KX) 

RMIN  =  0.  0 
R  M iX=  0 . 0 
DO  100  1=1, KX 

IPU(I)  .GT.  F  RAX)  RKAX=  X  (I) 

IF  * X  ( I)  .LT.  FEIN)  RKIN  =  X  (I) 

100  CONTINUE 
RETURN 
ENJ 
C 
C 

SUBROUTINE  PLX 8 AT  (Y 1 , Y 2 ,K1  Y  1 ,  IC2Y  1 ,  K  1Y2,K 2Y2,  X V  ECT  1,1  VECT 1 , 
IYV-CT2,  DT,ITIT1,JLAB1) 

DIMENSION  ITIT1  (10) 

DIN tN SION  JLABI(IC) 

DIMENSION  Y1  (K1Y1  ,K2Y  1) 

D1 JENSI ON  Y2(KlY2,K2Y2) 

DIMENSION  X VECT  1  ( K 2 Y 2 ) 

DIMENSION  YVBCT1  (K2Y2) 

DIJENSION  YVECT2  (K2Y2) 

DATA  ICBARA,1CHAFE/,A',(E*/ 

DO  20  I=1,K2Y2 
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X  V  SC  T 1  (I)  =DT*I 
20  CONTINUE 

DO  200  J*  1,  K IT 2 
DO  50  I=1,K2Y2 
I V  ECT  1  (I)  =Y  1  (J,I) 

Y VSCT 2(1)  *Y2  (J,I) 

50  CONTINUE 

CALL  VHNMX(YVECT1,K2Y2,Y1MIN,Y1HAX) 

CALL  VHNMX  (YVECI2  #K2Y2f  I2HIN,Y2MAX) 

IP ( Y 1  HI N  .LT.  Y2MIN)  YHIN*Y1MIN 
IP(Y2HIN  .LT.  Y  1  MIN) Y  HI  N*Y  2HI N 
I F  (  X  2  HI  N  .EC.  Y1MIN)  YHIN=Y2HIN 
IF ( Y 1 M AX  . GT.  Y2MAX) YHAX=Y 1  MAX 
IF (Y2HAX  .GI.  Y1HAX)  YHAX=T2HAX 
IF  ( Y2M A  X  .£0.  X 1 M AX) YHAX-Y2MAX 
CALL  WPLOT2  (X  VECT 1  (K2Y2),XVECT1  (1)  ,YMAX,YMIN) 

CA-L  MPLOT3  (IC  HA  BA  ,  XVE3 1 1 ,  YV  ECT  1  ,K2Y  2) 

CA-L  NPL0T3  (ICH  AN  E ,  X  V  E  Z II  ,  Y  V  ECT2  ,  K2  Y  2) 

WRITE  (6,90)  ITIT1 
90  POiMA  T  ( '  1'  ,40X,  10A4) 

C  A..L  WPLOT4  (40, JLAB1) 

200  CONTINUE 
fl  E  i"J  P  N 
ENj 

c 

c 

S  UoPOUTIN E  PLXSEN (DX H A T , K IX , K2X , X V ECT 1 , Y VECT 1 , DT, ITIT2 , JLA £2) 
DIMENSION  DXH AT  (K  1  X, K2  X) 

DIMENSION  XV  ECT 1  (K2X) 

DIMENSION  YV2CT1  ( K 2 X ) 

DIMENSION  I  TIT  2(10) 

DIMENSION  JLA B  2  ( 1  0) 

DO  20  1*1, K2X 
XVaCTI  (I) =DT*I 
20  CONTINUE 

DO  100  J=1,K1X 
DO  50  1*1, K2X 
YVaCTI  (I)  =DXH  AT  (J,I) 

50  CONTINUE 

WRITE  (6,60)  ITIT 2 
60  FOEMATC  1 ',40X,  10A4) 

CAjlL  WPLOT1  (XVECT1  ,YVECT1,K2X,  40,  JLAB2) 

100  CONTINUE 
RZIURN 
£Nu 
C 
C 

SUBROUTINE  CPICT (X, K 1 X , K2X, XVECT 1 , 1 V ECT  1 , THIN , IK  AX, DT, JLABEL, 
SITaTLE) 

DIMENSION  X  (K 1 X , R2X) 

DIMENSION  Y VECT  1  (K2X) 

DIMENSION  XVECT  1  (K2X) 

DIMENSION  JIABEL  (1  0) 

DIMENSION  I  TITLE  (10) 

DATA  ICHAR/'*'/ 

DO  20  I* 1 , K2X 
X  VaCT 1  (I) =DT*I 
20  CONTINUE 

DO  200  J=1,K1X 
DO  40  1*1, K2X 


f 


Ill 


YV2CT1  (I)  -I  (J,I) 

#0  CONTINUE 

WRITE (6,1 00) IT1TLE 
100  FORMAT {• 1 ##40X»  10A4) 

CAi.L  WPLOT2  (XV  ECT  1  (K2X)  ,XV2CT1(1)  .YKAI.YMIN) 
CALL  WPLOT3  (I  CH  A  R ,  XVLCT 1  #YVECT1,K2X) 

CALL  WPLOT4 (40, JLABEL) 

200  CONTINUE 
RETURN 
ENJ 
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Appendix  F 

LIST  OF  APL  FUNCTIONS  USED  TO  GENERATE  WHEELSET  SIMULATED  DATA 
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I 

I 

I 

I 

-|  *«*•*-» 
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C13 

C23 

C33 

C4] 

C5D 

C63 

C  73 

C83 

C93 

CIO] 

Cll] 

C12] 

C 13] 

C14] 

C 15] 

C  16] 

Cl  7] 

L 183 

C 1 9] 

C20] 

C21  ] 

C22] 

C23] 

C24] 


Cll 

C2] 

C3] 

C4] 

C5] 

C6] 

C7] 

C83 

C9] 

CIO] 


Cl] 
C2] 
L  3] 
C4] 


V  WHLSIM 
flt-SYSMTK 
Bf-XHF-MTX 

G<-  3  3  fO 
*«-  3  1  F  0 

H«-  3  3  f  100010001 

•RANDOM  INPUT  STATISTICS' 

Qf-  '  ENTER  SUM  LENGTH  MEAN  MNSQVALUE • 
STATRIf Q 

WfWHTHOISE  STATRI 

•DETERMINISTIC  INPUT  STATISTICS' 

Qf ' ENTER  SUM  LENGTH  MEAN  MNSOVALUE ' 
STATDIf D 

WT f WHTNO I SE  STATDI 

Uf-  (  STATRI  t  2]  ,  1  ,  DfW  +  WT 

D  IN 

;<MEA5>(- 

•MEASUREMENT  NOISE  STATISTICS* 

Qf  '  ENTER  SUM  LENGTH  MEAN  MHSCIVflLUE  ' 
STATMNf 0 

VlfWHTNOISE  STATMN 
V2fWHTNOISE  STATMN 
V3fWHTNOISE  STATMN 

v<-(3,statmnc2]  )f  vi  ,V2,  V3 

2f(H+, XNMEAS)+V 

V 


V  PfSYSMTX 

F1  1«-  ~1  X  (  (  (  2x  THETA  XCOEE22  )  +  <  *0X  VEL  )  )  -f  M  AS S  X  VEL  ) 
r12<-2xcoEF22x  theta  --mass 
p- 1  3<--  1  X  K  T-f  MASS 
F21«-0 

F22<-“1  XKSI  XVEL  f^XTHETAxCOEFUxl-Xl- 

F^SR-!  XALFHAxVEL-j-LxPO 

F31<-1 

F32«-0 

F33<-0 

F<  3  3  /‘Fll,F12»F13iF21»F22»F23»F31iF32»F33 
v 


V  GfINPMTX 

GlU-E-OfMASS 
G21  «-0 
G31<--1 

Cf  3  1  fGll »G21 tG31 
v 


I 

I 

V  YfWHTMO I  SE  K  j  A { 6 { MM  f MS 

I  Cl]  A«-?(xci]x;:C2])f  100000 

C2]  A«-<xC13»*C23>f  « 

C33  *«-<l»*C23)f+/« 

IC43  MNH+/t>T(x/fS) 

C53  »«-C1x(mm-xc3])  )+6 

[^]  MSfMNSOVflL  D 
.  C7]  *<-<  (5<C43^MS>*0»5)XB 

I  C83  i>* 

♦ 

?  Yf MN5GVAL  X 

Cl  3  »><+/(«+.  X<SK)  )r(X/(f!!)  ) 


y  y«-cov  ::;mh 

C 1  3  mn«-  ( «f  x )  f  ( +/:: )  *  “  1  tr  '< 

C23  MN«-$MN 

C  3  3  >•<-<<  «-MN )  + .  x at  ( K-MM ) )  .f- 1 1 f : : 


V  DYN 

Cl  3  TH-lfDTxiM 
[23  INITIALISE 
[33  F-HIfl'T  STM  A 
C4]  DM  AT  F:  I X 

C53  M«-0 

C63  lo:h<-m+i 

C73  >:*Ci  r  <M+l)3<-<FHi+,xx>:ci  iM])+DD4  .  xuuc;  ;m] 

C83  »-»Xf }  (M+i)n<-(c+,xx::Cf  f  <M+D3)  +  f+.xuucf  ;  (M+i)] 

C93  +(M<H-1)/L0 

CIO]  5 :«-<NN,N)f  ,XX 
Cl 1 3  Yf(LUrM)f,YY 

V 


V  INITIALIZE 

Cl  3  ))>u-(r  ,«>*0.5 

C2D  MMfrffcm 
C33  L-Lf-ffCc;  13 
C43  «««■  (MW»l»N)fO 

C53  Y»<-<LL,l,N)fO 
C63 

C73  ‘ENTER  K(0)‘ 

C83  ;:::C  51?  1 3<-Q 

C93  rrc  ;  l ;  i  n<-c+ .  xxx[  f  1  i  1 3 
v 
♦ 


V  F'HIfDT  STM  AJRJNJFACTORJ  INDEX 

Cl  3  <f  »A)*0.5 

C.23  i«-( \R)  o  ,  =  \R 

C33  ««-  0 

C43  FHifi 
[.53  FACT ORf  I 
C  A3  l-OJNfn+i 

C73  F ACTORf-F  ACTORf  ,  xAxl'TfN 
C83  F'H  I  i-  F'HI+F  ACTOR 

C93  INDEX*.  <+/ [  ,  FACTOR  )4-R*2 
L  103  •■»(  INDEX>EPS)/L0 

V 

♦ 


9  DMA  TP:  X  >1 1  H 
£13  H4-DT  +  4 

£23  dd<-  (N-f  3 )  *  ( <0  *™  A)+(4xm  stm  A)+(2x(2xh>  stm  A)4-(4x(3xM)  stm  a>  +  (DT 


Appendix  G 

LIST  OF  THE  PROGRAM  USED  TO  OPERATE  THE  A/D  CONVERTER 
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Appendix  H 

PLOTS  OF  THE  LIKELIHOOD  FUNCTION,  OBSERVATION  TERM  AND 
BIAS  TERM  FOR  TEST  CASES  1  THROUGH  9 


0.500  0.690  0.880  1.070  1.260  1.650 

Figure  H.l  Plot  of  the  likelihood  function  for  test  case  1. 


Fifnire  H.<?  riot,  of  the  '  rv.if.,  n  -rm  for  tent,  cane 


Figure  Hjt  Plot  of  the  likelihood  function  for  test  case 


Figure  H.8  Plot  of  the  observation  term  for  test  case 


Fif^ire  H.9  Plot  of  the  bias  term  for  ter.t  car.e 


LIKELIHOOD  EDUCTION 


Figure  H.10  Plot  of  the  likelihood  function  for  test  case 


Figure  H.ll  Plot  of  the  observation  term  for  test  ease  *>. 


Figure  H.13  Plot  of  the  likelihood  function  for  test  case 


Figure  H.lU  Plot  of  the  observation  term  for  test  ease 


Figure  H.15  Plot  of  the  bias  term  for  nest  ease 
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LIKELIHOOD  FUNCTION 
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Figure  H.19  Plot  of  the  likelih  od  function  for  tent  cane 


HKROR  TERR  OF  L 


Figure  H .20  Plot  of  the  observation  term  for  test  case 


Figure  4.21  Plot  of  the  bias  term  for  test  case 


Figure  H.22  Plot  of  the  likelihood  function  for  test  case 


Figure  H.23  Plot  of  the  observation  term  for  test  case 


LI Kl LI  HOOD  POHCTXO* 


Figure  H.25  Plot  of  the  likelihood  function  for  tent  case 


I 
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Figure  H.26  Plot  of  the  observation  term  for  test  case 


Figure  H.27  Plot  of  the  bias  term  for  tost  case 


